diff -urN ORG/gdmp/@mkall gdmp/@mkall
--- ORG/gdmp/@mkall	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/@mkall	2020-11-13 12:00:00.000000000 +0900
@@ -11,6 +11,7 @@
 fort  -Wall gadd.f90     -L/home/SHARE/lib -lgm -o ../bin/gadd
 fort  -Wall gsub.f90     -L/home/SHARE/lib -lgm -o ../bin/gsub
 fort  -Wall gtrim.f90    -L/home/SHARE/lib -lgm -o ../bin/gtrim
+fort  -Wall gflat.f90    -L/home/SHARE/lib -lgm -o ../bin/gflat
 fort  -Wall gtopo.f90    -L/home/SHARE/lib -lgm -o ../bin/gtopo
 fort  -Wall govlay.f90   -L/home/SHARE/lib -lgm -o ../bin/govlay
 fort  -Wall gojoin.f90   -L/home/SHARE/lib -lgm -o ../bin/gojoin
@@ -25,6 +26,7 @@
 fort  -Wall shade.f90    -L/home/SHARE/lib -lgm -o ../bin/shade
 fort  -Wall plmapl.f90   -L/home/SHARE/lib -lgm -o ../bin/plmapl
 fort  -Wall plmapcl.f90  -L/home/SHARE/lib -lgm -o ../bin/plmapcl
+fort  -Wall plxsec.f90   -L/home/SHARE/lib -lgm -o ../bin/plxsec
 fort  -Wall rearx.f90    -L/home/SHARE/lib -lgm -o ../bin/rearx
 fort  -Wall xplmap.f90   xplgobj.c  -L/home/SHARE/lib -lgm -o ../bin/xplmap
 fort  -Wall xplmapc.f90  xplgobj.c  -L/home/SHARE/lib -lgm -o ../bin/xplmapc
diff -urN ORG/gdmp/gflat.f90 gdmp/gflat.f90
--- ORG/gdmp/gflat.f90	1970-01-01 09:00:00.000000000 +0900
+++ gdmp/gflat.f90	2020-11-13 12:00:00.000000000 +0900
@@ -0,0 +1,106 @@
+!--------( GFLAT : gen. new (flat) grid-data of constant value )--------
+  implicit none
+  integer :: lwkdir, lrtrim
+  character(80) :: wdr, flog, fnam
+  character(20) :: sdtm;  character(8) :: narea, salt
+  integer :: ldr, i, j, i1, i2
+  integer :: ncn, nccn, iorgn, korgn, ispan, ispbn
+  integer :: ixl, iyl, mszxn, mszyn, mxn, myn, ixh, iyh
+  real :: xl, yl, xh, yh, val, vnul, alt
+
+  call premsg('GFLAT : gen. new (flat) grid-data of constant value')
+  call opnpin()
+  ldr = lwkdir(50,wdr) + 1
+  flog = wdr;  fnam = wdr
+  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
+  call premsg('')
+  call gparma('Enter Name-label for new grid ==> ', 8, narea)
+  call gparmi('Select Coordinate Number (+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 gparmi('  Central Longitude ?  (Deg.) : ', i1)
+    korgn = i1*60
+  else
+    if (ncn == 72) call abendm('nc=72 not supported')
+    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 premsg('Specify gridding parameters')
+  call premsg('  Southwest corner Coord. in km ?')
+  call gparmf('    Northing : ', xl)
+  ixl = nint(xl*1000.)
+  call gparmf('    Easting  : ', yl)
+  iyl = nint(yl*1000.)
+  call premsg('  Mesh size and numbers to N ?')
+  call gparmi('    Unit Mesh size in m :  ', mszxn)
+  if (mszxn <= 0) call abendm('invalid value')
+  call gparmi('    Number of Mesh to N :  ', mxn)
+  ixh = ixl + (mxn-1)*mszxn
+  xh = float(ixh) / 1000.
+  call premsg('  Mesh size and numbers to E ?')
+  call gparmi('    Unit Mesh size in m :  ', mszyn)
+  if (mszyn <= 0) call abendm('invalid value')
+  call gparmi('    Number of Mesh to E :  ', myn)
+  iyh = iyl + (myn-1)*mszyn
+  yh = float(iyh) / 1000.
+  call gparmf('  (Constant) Mesh value ? ==> ', val)
+  if ((val > 9999.) .or. (val < -9999.))   &
+    call abendm('MeshVal: Out of valid range')
+  vnul = 9999.9
+  call premsg('  Enter const. Obs. Altitude (> 0.)')
+  call gparmf('    (if 0 or minus, Undef. assumed) ==> ', alt)
+  if (alt > 9999.) call abendm('Alt: Out of valid range')
+  if (alt > 0.) then
+    write(salt,'(1x,f7.0)') alt
+  else
+    alt = -1.;  salt = '  Undef.'
+  endif
+  call gparma('Enter new grid output filename ==> ', 81-ldr, fnam(ldr:80))
+  call clspin()
+  open(2,file=fnam,status='new')
+  write(6,'(a)') '----------------------------------------'
+  write(6,'(a)') 'GFLAT : gen. new (flat) grid-data of constant value'
+  write(6,'(a,a)') ' Output new file : ', fnam(1:lrtrim(fnam))
+  write(6,'(a12,a8,6x,a6,i4,4x,4i8)')   &
+    'Areaname :  ', narea, 'Coord.', nccn, iorgn,korgn, ispan,ispbn
+  write(6,'(4x,a)')   &
+    '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
+  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
+    'gen.new grid', xl, yl, mszxn, mszyn, mxn, myn, salt
+
+  write(2,'(a)')  '# prog.GFLAT output'
+  write(2,'(a8,4x,i4,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn
+  write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixl,iyl, mszxn,mszyn, mxn,myn, vnul, alt
+  do j=1,myn
+    write(2,'((f7.1,9(1x,f7.1)))') (val,i=1,mxn)
+  enddo
+  close(2)
+
+  if (flog(ldr:ldr) /= ' ') then
+    open(99,file=flog,access='append')
+    write(99,'(//a,$)') 'Gflat :'
+    write(99,'(a)')   ' gen. new (flat) grid-data of constant value'
+    write(99,'(a,a)') ' Output new file : ', fnam(1:lrtrim(fnam))
+    write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
+      'Areaname :  ', narea, 'Coord.', nccn, iorgn,korgn, ispan,ispbn
+    write(99,'(4x,a)')   &
+      '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
+    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
+      'gen.new grid', xl, yl, mszxn, mszyn, mxn, myn, salt
+    call strdtm(sdtm)
+    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
+    close(99)
+  endif
+  stop
+  end
diff -urN ORG/gdmp/plmap.f90 gdmp/plmap.f90
--- ORG/gdmp/plmap.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/plmap.f90	2019-12-17 10:44:00.000000000 +0900
@@ -4,9 +4,9 @@
   real,allocatable :: g(:)
   character(80) :: wdr, flog, fnin, fnam, buf
   character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
-  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
+  character(6) :: ss;  character(3) :: spl = '   ';  character(1) :: yn
   integer :: ldr, mc, i, j, n, k, IER
-  integer :: istep, istp, lcap, kpl, kscl, klll, kshl, km, iinc, kinc
+  integer :: istep, istp, lcap, kpl, kscl, klll, kshl, km, ls, iinc, kinc
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: isf, itf, ksf, ktf, idf, kdf
@@ -125,17 +125,55 @@
       fscl = fscl * 100.
       ccm = fscl/100000.
       kscl = 0;  sclr = 0.;  sclu = 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
@@ -204,9 +242,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)
diff -urN ORG/gdmp/plmapc.f90 gdmp/plmapc.f90
--- ORG/gdmp/plmapc.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/plmapc.f90	2019-12-17 10:44:00.000000000 +0900
@@ -4,9 +4,9 @@
   real,allocatable :: g(:)
   character(80) :: wdr, flog, fnin, fnam, buf
   character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
-  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
+  character(6) :: ss;  character(3) :: spl = '   ';  character(1) :: yn
   integer :: ldr, mc, i, j, n, k, IER
-  integer :: istep, istp, kpl, kscl, klll, kshl, km, iinc, kinc
+  integer :: istep, istp, kpl, kscl, klll, kshl, km, ls, iinc, kinc
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: isf, itf, ksf, ktf, idf, kdf, mval
@@ -125,17 +125,55 @@
       fscl = fscl * 100.
       ccm = fscl/100000.
       kscl = 0;  sclr = 0.;  sclu = 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
@@ -206,9 +244,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)
diff -urN ORG/gdmp/plmapcl.f90 gdmp/plmapcl.f90
--- ORG/gdmp/plmapcl.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/plmapcl.f90	2021-08-11 20:00:00.000000000 +0900
@@ -4,14 +4,15 @@
   real,allocatable :: g(:)
   character(80) :: wdr, flog, fnin, fnln, fnam, buf
   character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
-  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
+  character(6) :: ss;  character(3) :: spl = '   ';  character(1) :: yn
   integer :: ldr, mc, i, j, n, k, IER
-  integer :: istep, istp, kpl, kscl, klll, kshl, km, iinc, kinc
+  integer :: istep, istp, kpl, kscl, klll, kshl, km, ls, iinc, kinc
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: isf, itf, ksf, ktf, idf, kdf, kl, mval, np
-  real :: vnul, xs, ys, gsum, gmin, gmax, gdev, dg, alat, alon, vhi, vlo
-  real :: fscl, high, wide, ccm, sclr, sclu, dscl, xo, yo, xt, yt, x, y
+  real(8) :: dlat, dlon, dx, dy
+  real :: vnul, xs, ys, gsum, gmin, gmax, gdev, dg, vhi, vlo
+  real :: fscl, high, wide, ccm, sclr, sclu, dscl, xo, yo, xt, yt
   real :: fi0, fi1, fi2, fi3, fi4, fi5, fk0, fk1, fk2, fk3, fk4, fk5
   external cviken
 
@@ -131,17 +132,55 @@
       fscl = fscl * 100.
       ccm = fscl/100000.
       kscl = 0;  sclr = 0.;  sclu = 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
@@ -213,9 +252,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)
@@ -239,21 +277,21 @@
         do i=isf/iinc,itf/iinc
           idf = i * iinc
           if (idf < isf.or.idf > itf) cycle
-          call cviken(float(idf), float(ksf), y, x)
-          call plot((y-ys)/ccm, (x-xs)/ccm, 3)
+          call cvdiken(dble(idf), dble(ksf), dy, dx)
+          call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 3)
           do kdf=ksf+1,ktf
-            call cviken(float(idf), float(kdf), y, x)
-            call plot((y-ys)/ccm, (x-xs)/ccm, 2)
+            call cvdiken(dble(idf), dble(kdf), dy, dx)
+            call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 2)
           enddo
         enddo
         do k=ksf/kinc,ktf/kinc
           kdf = k * kinc
           if (kdf < ksf.or.kdf > ktf) cycle
-          call cviken(float(isf), float(kdf), y, x)
-          call plot((y-ys)/ccm, (x-xs)/ccm, 3)
+          call cvdiken(dble(isf), dble(kdf), dy, dx)
+          call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 3)
           do idf=isf+1,itf
-            call cviken(float(idf), float(kdf), y, x)
-            call plot((y-ys)/ccm, (x-xs)/ccm, 2)
+            call cvdiken(dble(idf), dble(kdf), dy, dx)
+            call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 2)
           enddo
         enddo
       endif
@@ -272,9 +310,9 @@
           if (i /= 0) buf(i:i) = ' '
           i = index(buf,'E')
           if (i /= 0) buf(i:i) = ' '
-          read(buf,*) alat, alon
-          call cviken(alat, alon, y, x)
-          call plot((y-ys)/ccm, (x-xs)/ccm, 3)
+          read(buf,*) dlat, dlon
+          call cvdiken(dlat, dlon, dy, dx)
+          call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 3)
           kl = 1
         endif
         read(11,'(a)',iostat=IER) buf
@@ -286,9 +324,9 @@
         if (i /= 0) buf(i:i) = ' '
         i = index(buf,'E')
         if (i /= 0) buf(i:i) = ' '
-        read(buf,*) alat, alon
-        call cviken(alat, alon, y, x)
-        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
+        read(buf,*) dlat, dlon
+        call cvdiken(dlat, dlon, dy, dx)
+        call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 2)
       enddo
       rewind(11)
       call newpen(1)
diff -urN ORG/gdmp/plmapcs.f90 gdmp/plmapcs.f90
--- ORG/gdmp/plmapcs.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/plmapcs.f90	2019-12-17 10:44:00.000000000 +0900
@@ -5,9 +5,9 @@
   real,parameter :: RAD = (180./3.141593)
   character(80) :: wdr, flog, fnin, fnam, buf
   character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
-  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
+  character(6) :: ss;  character(3) :: spl = '   ';  character(1) :: yn
   integer :: ldr, mc, i, j, n, k, IER
-  integer :: istep, istp, kpl, kscl, klll, kshl, km, iinc, kinc
+  integer :: istep, istp, kpl, kscl, klll, kshl, km, ls, iinc, kinc
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: isf, itf, ksf, ktf, idf, kdf, mval
@@ -137,17 +137,55 @@
       fscl = fscl * 100.
       ccm = fscl/100000.
       kscl = 0;  sclr = 0.;  sclu = 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
@@ -264,9 +302,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)
diff -urN ORG/gdmp/plmapg.f90 gdmp/plmapg.f90
--- ORG/gdmp/plmapg.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/plmapg.f90	2019-12-17 10:44:00.000000000 +0900
@@ -4,9 +4,9 @@
   real,allocatable :: g(:)
   character(80) :: wdr, flog, fnin, fnam, buf
   character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
-  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
+  character(6) :: ss;  character(3) :: spl = '   ';  character(1) :: yn
   integer :: ldr, mc, i, j, n, k, IER
-  integer :: istep, istp, kpl, kscl, klll, kshl, km, iinc, kinc
+  integer :: istep, istp, kpl, kscl, klll, kshl, km, ls, iinc, kinc
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: isf, itf, ksf, ktf, idf, kdf, mval
@@ -125,17 +125,55 @@
       fscl = fscl * 100.
       ccm = fscl/100000.
       kscl = 0;  sclr = 0.;  sclu = 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
@@ -204,9 +242,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)
diff -urN ORG/gdmp/plmapl.f90 gdmp/plmapl.f90
--- ORG/gdmp/plmapl.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/plmapl.f90	2021-08-11 20:00:00.000000000 +0900
@@ -4,14 +4,15 @@
   real,allocatable :: g(:)
   character(80) :: wdr, flog, fnin, fnln, fnam, buf
   character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
-  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
+  character(6) :: ss;  character(3) :: spl = '   ';  character(1) :: yn
   integer :: ldr, mc, i, j, n, k, IER
-  integer :: istep, istp, lcap, kpl, kscl, klll, kshl, km, iinc, kinc
+  integer :: istep, istp, lcap, kpl, kscl, klll, kshl, km, ls, iinc, kinc
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: isf, itf, ksf, ktf, idf, kdf, kl, np
-  real :: vnul, xs, ys, gsum, gmin, gmax, gdev, dg, alat, alon
-  real :: fscl, high, wide, ccm, sclr, sclu, dscl, xo, yo, xt, yt, x, y
+  real(8) :: dlat, dlon, dx, dy
+  real :: vnul, xs, ys, gsum, gmin, gmax, gdev, dg
+  real :: fscl, high, wide, ccm, sclr, sclu, dscl, xo, yo, xt, yt
   real :: fi0, fi1, fi2, fi3, fi4, fi5, fk0, fk1, fk2, fk3, fk4, fk5
   external cviken
 
@@ -131,17 +132,55 @@
       fscl = fscl * 100.
       ccm = fscl/100000.
       kscl = 0;  sclr = 0.;  sclu = 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
@@ -210,9 +249,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)
@@ -236,21 +274,21 @@
         do i=isf/iinc,itf/iinc
           idf = i * iinc
           if ((idf < isf) .or. (idf > itf)) cycle
-          call cviken(float(idf), float(ksf), y, x)
-          call plot((y-ys)/ccm, (x-xs)/ccm, 3)
+          call cvdiken(dble(idf), dble(ksf), dy, dx)
+          call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 3)
           do kdf=ksf+1,ktf
-            call cviken(float(idf), float(kdf), y, x)
-            call plot((y-ys)/ccm, (x-xs)/ccm, 2)
+            call cvdiken(dble(idf), dble(kdf), dy, dx)
+            call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 2)
           enddo
         enddo
         do k=ksf/kinc,ktf/kinc
           kdf = k * kinc
           if ((kdf < ksf) .or. (kdf > ktf)) cycle
-          call cviken(float(isf), float(kdf), y, x)
-          call plot((y-ys)/ccm, (x-xs)/ccm, 3)
+          call cviken(dble(isf), dble(kdf), dy, dx)
+          call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 3)
           do idf=isf+1,itf
-            call cviken(float(idf), float(kdf), y, x)
-            call plot((y-ys)/ccm, (x-xs)/ccm, 2)
+            call cviken(dble(idf), dble(kdf), dy, dx)
+            call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 2)
           enddo
         enddo
       endif
@@ -270,9 +308,9 @@
           if (i /= 0) buf(i:i) = ' '
           i = index(buf,'E')
           if (i /= 0) buf(i:i) = ' '
-          read(buf,*) alat, alon
-          call cviken(alat, alon, y, x)
-          call plot((y-ys)/ccm, (x-xs)/ccm, 3)
+          read(buf,*) dlat, dlon
+          call cvdiken(dlat, dlon, dy, dx)
+          call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 3)
           kl = 1
         endif
         read(11,'(a)',iostat=IER) buf
@@ -284,9 +322,9 @@
         if (i /= 0) buf(i:i) = ' '
         i = index(buf,'E')
         if (i /= 0) buf(i:i) = ' '
-        read(buf,*) alat, alon
-        call cviken(alat, alon, y, x)
-        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
+        read(buf,*) dlat, dlon
+        call cvdiken(dlat, dlon, dy, dx)
+        call plot(sngl(dy-ys)/ccm, sngl(dx-xs)/ccm, 2)
       enddo
       rewind(11)
       call newpen(1)
diff -urN ORG/gdmp/plmaps.f90 gdmp/plmaps.f90
--- ORG/gdmp/plmaps.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/plmaps.f90	2019-12-17 10:44:00.000000000 +0900
@@ -5,9 +5,9 @@
   real,parameter :: RAD = (180./3.141593)
   character(80) :: wdr, flog, fnin, fnam, buf
   character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
-  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
+  character(6) :: ss;  character(3) :: spl = '   ';  character(1) :: yn
   integer :: ldr, mc, i, j, n, k, IER
-  integer :: istep, istp, lcap, kpl, kscl, klll, kshl, km, iinc, kinc
+  integer :: istep, istp, lcap, kpl, kscl, klll, kshl, km, ls, iinc, kinc
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: isf, itf, ksf, ktf, idf, kdf
@@ -137,17 +137,55 @@
       fscl = fscl * 100.
       ccm = fscl/100000.
       kscl = 0;  sclr = 0.;  sclu = 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
@@ -266,9 +304,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)
diff -urN ORG/gdmp/plxsec.f90 gdmp/plxsec.f90
--- ORG/gdmp/plxsec.f90	1970-01-01 09:00:00.000000000 +0900
+++ gdmp/plxsec.f90	2020-11-21 12:30:00.000000000 +0900
@@ -0,0 +1,1034 @@
+!--------( PLXSEC : Plot Arbitrary cross-section of anomaly map )--------
+  implicit none
+  integer :: lwkdir, lrtrim
+  real,allocatable :: g(:,:), h(:,:), gs(:), hs(:)
+  real :: pam(20), pth(20), pph(20), pv(7,20)
+  real :: ds1(20), ds2(20), hi1(20), hi2(20)
+  integer :: mkp(22), jj(20)
+  character(80) :: wdr, flog, fnin, fnam, buf
+  character(72) :: out;  character(20) :: sdtm;
+  character(8) :: area, area1, unit;
+  character(6) :: ss;    character(3) :: shv;    character(1) :: yn
+  character(8) :: salt;  character(25) :: sloc;  character(10) :: slay
+  integer :: mk, ldr, mc, i, j, n, k, kloc, ndp, IOS
+  integer :: istep, istp, khv, kscl, klll, kshl, km, ls, iinc, kinc
+  integer :: nc, ncc, iorg, korg, ispa, ispb
+  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
+  integer :: ixs, iys, mszx, mszy, nmx, nmy
+  integer :: mdl=0, nmd=0, kmd, ks, isf, itf, ksf, ktf, idf, kdf, mval
+  integer :: ix, iy, ix1, iy1, mszx1, mszy1, nmx1, nmy1
+  real :: am, the, phi, v1, v2, v3, v4, v5, v6, v7
+  real :: alt, vnul, hnul, xs, ys, dx, dy, vlo, vhi, fx, fy
+  real :: gsum, gmin, gmax, gdev, dg, g1, g2, gb, gt, gr
+  real :: hsum, hmin, hmax, hdev, dh, h1, h2, hb, ht, hr
+  real :: gssum, gsmin, gsmax, hssum, hsmin, hsmax
+  real :: fscl, ccm, sclr, sclu, dscl, xo, yo, xt, yt, x, y
+  real :: fi0, fi1, fi2, fi3, fi4, fi5, fk0, fk1, fk2, fk3, fk4, fk5
+  real :: high, him, hig, hih, wide, wim, wig
+  real :: p1, p2, p3, q1, q2, q3, smx, smy, xinc, yinc
+  real :: xp, yp, xpa, ypa, xpb, ypb, xst, yst, xed, yed
+  real :: xna, yea, xnb, yeb, dmy, dip, dec
+  real :: deep, xc, xr, yc, yr, za, zr, ds, pr
+  external cviken
+
+  call premsg('PLXSEC : Plot arbitrary cross-section of anomaly map')
+  call opnpin()
+  ldr = lwkdir(50,wdr) + 1
+  flog = wdr;  fnin = wdr;  fnam = wdr
+  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
+  call gparma('Enter Input data filename ==> ', 81-ldr, fnin(ldr:80))
+  open(10,file=fnin,status='old')
+  call gparma('Enter Output PS filename  ==> ', 81-ldr, fnam(ldr:80))
+  if (fnam(ldr:ldr) == ' ') fnam(1:1) = char(0)
+  if (fnam(80:80) /= ' ') call abendm('too long filename')
+
+  read(10,'(a)') buf
+  do while (buf(1:1) == '#')
+    if (mdl == 0) then
+      if (buf(1:20) == '### CALMAS model ###') mdl = 1
+    else if (mdl == 1) then
+      if ((buf(1:12)=='##Block  1  ') .or. (buf(1:12)=='##hRect  2  ') .or.   &
+          (buf(1:12)=='##vLine  3  ') .or. (buf(1:12)=='##Point  4  ')) then
+        read(buf(8:80),*,IOSTAT=IOS) mk, am, the, phi
+        mdl = 2;  if (IOS /= 0) mdl = -1
+      else if (buf(1:12) == '##AmbMF  0  ') then
+        read(buf(12:80),*,IOSTAT=IOS) dmy, dip, dec
+        mdl = 3;  if (IOS /= 0) mdl = -1
+      else
+        mdl = -1
+      endif
+    else if (mdl == 2) then
+      if (buf(1:4) == '#   ') then
+        read(buf(2:80),*,IOSTAT=IOS) v1,v2,v3,v4,v5,v6,v7
+        if (IOS == 0) then
+          if (nmd < 20) then
+            nmd = nmd + 1; mkp(nmd) = mk; pam(nmd) = am; pth(nmd) = the
+            pph(nmd) = phi; pv(1,nmd) = v1; pv(2,nmd) = v2; pv(3,nmd) = v3
+            pv(4,nmd) = v4; pv(5,nmd) = v5; pv(6,nmd) = v6; pv(7,nmd) = v7
+            mdl = 1
+          else
+            mdl = 4             !--- too many models
+          endif
+        else
+          mdl = -1
+        endif
+      endif
+    endif
+    read(10,'(a)') buf
+  enddo
+  if (mdl == 0) then
+    call premsg('--- CALMAS model inf. not found')
+  else if (mdl == 3) then
+    call premsg('--- CALMAS model inf. accepted')
+  else if (mdl == 4) then
+    call premsg('--- too many models, truncated to 20')
+  else
+    mdl = -1;  call premsg('--- invalid CALMAS model inf. neglected')
+  endif
+
+  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;  mc = 0
+    else
+      ncc = ncc - 200;  nc = ncc;  mc = 1
+    endif
+  else                                !-- v2018
+    read(buf,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
+    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
+    korg = (nc-30)*360 - 180
+    if (nc > 60) korg = 0
+    if (nc == 61) iorg = +5400
+    if (nc == 62) iorg = -5400
+  endif
+  read(10,*) ixs, iys, mszx, mszy, nmx, nmy, vnul, alt
+  xs = float(ixs)/1000.;  ys = float(iys)/1000.
+  write(salt,'(1x,f7.0)') alt
+  if (alt == 0.) salt = '    var.'
+  if (alt < 0.)  call abendm('Undefined Altitude')
+  if ((nmx <= 0) .or. (nmy <= 0)) call abendm('illegal array size')
+  allocate(g(nmy,nmx), h(nmy,nmx))
+  read(10,*) ((g(j,i),i=1,nmx),j=1,nmy)
+  n = 0;  k = 0
+  do i=1,nmx
+    do j=1,nmy
+      if (g(j,i) == vnul) cycle
+      gsum = 0.;  gmin = g(j,i);  gmax = g(j,i);  k = 1;  exit
+    enddo
+    if (k == 1) exit
+  enddo
+  if (k == 0) call abendm('no valid data')
+  do i=1,nmx; do j=1,nmy
+    if (g(j,i) == vnul) cycle
+    n = n + 1
+    gsum = gsum + g(j,i)
+    if (g(j,i) < gmin) gmin = g(j,i)
+    if (g(j,i) > gmax) gmax = g(j,i)
+  enddo; enddo
+  gsum = gsum / float(n)
+  n = 0;  gdev = 0.
+  do i=1,nmx; do j=1,nmy
+    if (g(j,i) == vnul) cycle
+    dg = g(j,i) - gsum
+    n = n + 1
+    gdev = gdev + dg*dg
+  enddo; enddo
+  gdev = sqrt(gdev / float(n))
+
+  if (alt == 0.) then
+    do
+      read(10,'(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
+    if ((area1 /= area) .or. (ncc1 /= ncc) .or. (iorg1 /= iorg)   &
+        .or. (korg1 /= korg) .or. (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
+      call abendm('Header1 inconsistent')
+    read(10,*) ix1, iy1, mszx1, mszy1, nmx1, nmy1, hnul
+    if ((ix1 /= ixs) .or. (iy1 /= iys) .or. (mszx1 /= mszx)   &
+        .or. (mszy1 /= mszy) .or. (nmx1 /= nmx) .or. (nmy1 /= nmy))   &
+      call abendm('Header2 inconsistent')
+    read(10,*) ((h(j,i),i=1,nmx),j=1,nmy)
+    n = 0;  k = 0
+    do i=1,nmx
+      do j=1,nmy
+        if (h(j,i) == hnul) cycle
+        hsum = 0.;  hmin = h(j,i);  hmax = h(j,i);  k = 1;  exit
+      enddo
+      if (k == 1) exit
+    enddo
+    if (k == 0) call abendm('no valid data')
+    do i=1,nmx; do j=1,nmy
+      if (h(j,i) == hnul) cycle
+      n = n + 1;  hsum = hsum + h(j,i)
+      if (h(j,i) < hmin) hmin = h(j,i)
+      if (h(j,i) > hmax) hmax = h(j,i)
+    enddo; enddo
+    hsum = hsum / float(n)
+    n = 0;  hdev = 0.
+    do i=1,nmx; do j=1,nmy
+      if (h(j,i) == hnul) cycle
+      dh = h(j,i) - hsum
+      n = n + 1
+      hdev = hdev + dh*dh
+    enddo; enddo
+    hdev = sqrt(hdev / float(n))
+  else
+    hnul = 9999.9
+    do i=1,nmx; do j=1,nmy
+      if (g(j,i) == vnul) then
+        h(j,i) = hnul
+      else
+        if (alt < 0.) then
+          h(j,i) = hnul
+        else
+          h(j,i) = alt
+        endif
+      endif
+    enddo; enddo
+    hmin = alt;  hmax = alt;  hsum = alt;  hdev = 0.
+  endif
+  close(10)
+  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
+    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
+  call premsg(out)
+  write(out,'(14x,2a10,2x,2a8,2a6,4x,a4)')   &
+    'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
+  call premsg(out)
+  write(out,'(16x,2f10.2,2i8,2i6,a8)') xs, ys, mszx, mszy, nmx, nmy, salt
+  call premsg(out)
+  call premsg('')
+  write(out,'(4x,3(a,f7.1),a,f6.1)') 'Data Range : [', gmin,   &
+    ' ,', gmax, ' ]   Mean :', gsum, '   St.Dev.:', gdev
+  call premsg(out)
+  if (alt == 0.) then
+    write(out,'(4x,3(a,f7.1),a,f6.1)') 'Alt. Range : [', hmin,   &
+      ' ,', hmax, ' ]   Mean :', hsum, '   St.Dev.:', hdev
+  else if (alt < 0.) then
+    write(out,'(4x,3a)') 'Obs. Altitude :  undefined'
+  else
+    write(out,'(4x,3a)') 'Obs. Altitude :  ', salt, ' m const.'
+  endif
+  call premsg(out)
+  call premsg('')
+
+  do i=1,nmd                   ! convert Coord. val to Dist. from SW-corner
+    pv(1,i) = pv(1,i) - xs;  pv(3,i) = pv(3,i) - ys
+  enddo
+
+  call premsg('[Section A-B definition]')
+  call premsg('  How to specify ends location')
+  call premsg('    0: Coordinates values   or')
+  call gparmi('    1: Distance from SW corner ==> ', kloc)
+  if (kloc == 0) then
+    sloc = 'Coordinate (km) value    '
+    call premsg('  Specify by Northing and Easting coordinates')
+    call prompt('    N and E coord. (in km) of (A) start point')
+    call gparmf2(' ==> ', xna, yea)
+    call prompt('    N and E coord. (in km) of (B) end point  ')
+    call gparmf2(' ==> ', xnb, yeb)
+    xpa = xna - xs;  ypa = yea - ys;  xpb = xnb - xs;  ypb = yeb - ys
+  else if (kloc == 1) then
+    sloc = 'Dist. (km) from SW corner'
+    call premsg('  Specify by relative distances from SW corner')
+    call prompt('    N and E distance (in km) of (A) start point')
+    call gparmf2(' ==> ', xpa, ypa)
+    call prompt('    N and E distance (in km) of (B) end point  ')
+    call gparmf2(' ==> ', xpb, ypb)
+    xna = xpa + xs;  yea = ypa + ys;  xnb = xpb + xs;  yeb = ypb + ys
+  else
+    call abendm('invalid selection')
+  endif
+  ds = sqrt((xpb-xpa)*(xpb-xpa) + (ypb-ypa)*(ypb-ypa))
+  if (ds < 0.01) call abendm('too short profile')
+  call premsg('  The section profile is approximated by a series of points')
+  call gparmi('    Enter no. of points consisting of the section ==> ', ndp)
+  call premsg('  Show source Location (if such data available) ?')
+  call gparma('        "y" for Yes, othewise No ==> ', 1, yn)
+  kmd = 0;  if ((yn == 'y') .or. (yn == 'Y')) kmd = 1
+  if (kmd == 1) then
+    if (mdl == 0) then
+      call premsg('--- CALMAS model inf. not found')
+    else if (mdl == 4) then
+      call premsg('--- too many models, truncated to 20')
+    else if (mdl == -1) then
+      call premsg('--- invalid CALMAS model inf. neglected')
+    endif
+  endif
+
+!-- Generate section data
+  allocate (gs(ndp), hs(ndp))
+  smx = float(mszx);  smy = float(mszy)
+  xp = xpa*1000.;  xinc = (xpb-xpa)*1000./(ndp-1)
+  yp = ypa*1000.;  yinc = (ypb-ypa)*1000./(ndp-1)
+  xst = xpa;  yst = ypa;  xed = xpb;  yed = ypb
+  do i=1,ndp
+    fx = xp/smx;  ix = ifix(fx);  dx = fx - float(ix)
+    if ((ix==(nmx-1)) .and. (dx==0.)) then
+      ix = nmx - 2; dx = 1.
+    endif
+    fy = yp/smy;  iy = ifix(fy);  dy = fy - float(iy)
+    if ((iy==(nmy-1)) .and. (dy==0.)) then
+      iy = nmy - 2; dy = 1.
+    endif
+    if ((ix<0) .or. (dx<0.) .or. (iy<0) .or.(dy<0.)) then
+      gs(i) = vnul;  hs(i) = hnul
+    else if ((ix>=(nmx-1)) .or. (iy>=(nmy-1))) then
+      gs(i) = vnul;  hs(i) = hnul
+    else if (g(iy+1,ix+1) == vnul) then
+      gs(i) = vnul;  hs(i) = hnul
+    else if (dy == 0.) then
+      if (dx == 0.) then
+        gs(i) = g(iy+1,ix+1);  hs(i) = h(iy+1,ix+1)
+      else if (g(iy+1,ix+2) == vnul) then
+        gs(i) = vnul;  hs(i) = hnul
+      else
+        gs(i) = g(iy+1,ix+1) + (g(iy+1,ix+2)-g(iy+1,ix+1))*dx
+        if ((h(iy+1,ix+1)==hnul) .or. (h(iy+1,ix+2)==hnul)) then
+          hs(i) = hnul
+        else
+          hs(i) = h(iy+1,ix+1) + (h(iy+1,ix+2)-h(iy+1,ix+1))*dx
+        endif
+      endif
+    else if (g(iy+1,ix+2) == vnul) then
+      gs(i) = vnul;  hs(i) = hnul
+    else
+      if (dx == 0.) then
+        gs(i) = g(iy+1,ix+1) + (g(iy+2,ix+1)-g(iy+1,ix+1))*dy
+        if ((h(iy+1,ix+1)==hnul) .or. (h(iy+2,ix+1)==hnul)) then
+          hs(i) = hnul
+        else
+          hs(i) = h(iy+1,ix+1) + (h(iy+2,ix+1)-h(iy+1,ix+1))*dy
+        endif
+      else if ((g(iy+2,ix+1)==vnul) .or. (g(iy+2,ix+2)==vnul)) then
+        gs(i) = vnul;  hs(i) = hnul
+      else
+        g1 = g(iy+1,ix+1) + (g(iy+1,ix+2)-g(iy+1,ix+1))*dx
+        g2 = g(iy+2,ix+1) + (g(iy+2,ix+2)-g(iy+2,ix+1))*dx
+        gs(i) = g1 + (g2-g1)*dy
+        if ((h(iy+1,ix+1)==hnul) .or. (h(iy+1,ix+2)==hnul) .or.   &
+            (h(iy+1,ix+1)==hnul) .or. (h(iy+2,ix+2)==hnul)) then
+          hs(i) = hnul
+        else
+          h1 = h(iy+1,ix+1) + (h(iy+1,ix+2)-h(iy+1,ix+1))*dx
+          h2 = h(iy+2,ix+1) + (h(iy+2,ix+2)-h(iy+2,ix+1))*dx
+          hs(i) = h1 + (h2-h1)*dy
+        endif
+      endif
+    endif
+    xp = xp + xinc;  yp = yp + yinc
+  enddo
+  n = 0
+  do i=1,ndp
+    if (gs(i) == vnul) cycle
+    gssum = gs(i);  gsmin = gs(i);  gsmax = gs(i);  n = 1;  exit
+  enddo
+  if (n == 0) call abendm('no valid data')
+  do j=i+1,ndp
+    if (gs(j) == vnul) cycle
+    n = n + 1
+    gssum = gssum + gs(j)
+    if (gs(j) < gsmin) gsmin = gs(j)
+    if (gs(j) > gsmax) gsmax = gs(j)
+  enddo
+  gssum = gssum / float(n)
+  n = 0
+  do i=1,ndp
+    if (hs(i) == hnul) cycle
+    hssum = hs(i);  hsmin = hs(i);  hsmax = hs(i);  n = 1;  exit
+  enddo
+  if (n == 0) call abendm('no valid data')
+  do j=i+1,ndp
+    if (hs(j) == hnul) cycle
+    n = n + 1
+    hssum = hssum + hs(j)
+    if (hs(j) < hsmin) hsmin = hs(j)
+    if (hs(j) > hsmax) hsmax = hs(j)
+  enddo
+  hssum = hssum / float(n)
+
+  call premsg('[Map-Section Layout / Paper orientation]')
+  call premsg('  Select Layout  1: horizontal  (landscape)')
+  call gparmi('            or   2: vertical    (portrait)  ==> ', khv)
+  if (khv == 1) then
+    shv = 'A4L';  slay = 'Horizontal'
+  else if (khv == 2) then
+    shv = 'A4P';  slay = '  Vertical'
+  else
+    call abendm('invalid value')
+  endif
+  call premsg('[Index map]')
+  if (khv == 1) then
+    call premsg('    (on the Left side aligned Middle)')
+    call gparmf('  Specify Width of drawing Index map (in cm) ==> ', wim)
+    fscl = float((nmy-1)*mszy) / wim
+    him = float((nmx-1)*mszx) / fscl
+    write(out,'(4x,a,f8.1,a1,f8.1,a1)')  &
+      'Index map size :', wim, 'W', him, 'H'
+  else
+    call premsg('    (on the Top side aligned Center)')
+    call gparmf('  Specify Height of drawing Index map (in cm) ==> ', him)
+    fscl = float((nmx-1)*mszx) / him
+    wim = float((nmy-1)*mszy) / fscl
+    write(out,'(4x,a,f8.1,a1,f8.1,a1)')  &
+      'Index map size :', him, 'H', wim, 'W'
+  endif
+  fscl = fscl * 100.
+  call premsg('  Specify contour interval (nT)')
+  call gparmi('    (1/2 of color-grading interval) ==> ', istep)
+  if (istep == 0) call abendm('invalid value')
+  istp = iabs(istep)
+  call gparmi('  Median value of color-grading ?  ==> ', mval)
+  vlo = float(mval - istp*20)
+  vhi = float(mval + istp*20)
+  ccm = fscl/100000.
+  kscl = 0;  sclr = 0.;  sclu = 0.;  klll = 0;  kshl = 0
+  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
+    kscl = 1
+    write(out,'(a,f4.1,a)') '      size :', dscl, 'cm Wide, 0.8cm High'
+    call premsg(out)
+    call prompt('    Lower-Left Pos. ? [right and up from SW corner]')
+    call gparmf2(' (in cm) ==> ', sclr, sclu)
+  endif
+  call prompt('  Draw Lat.Lon.Lines ?')
+  call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
+  if ((yn == 'y') .or. (yn == 'Y')) klll = 1
+  call prompt('  Draw Coastlines ?   ')
+  call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
+  if ((yn == 'y') .or. (yn == 'Y')) then
+    kshl = 2
+    call prompt('       Rivers  ?      ')
+    call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
+    if ((yn == 'y') .or. (yn == 'Y')) kshl = 12
+    call prompt('       Pref.Boundary ?')
+    call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
+    if ((yn == 'y') .or. (yn == 'Y')) then
+      kshl = kshl + 200
+      if (kshl == 202) kshl = 102
+    endif
+  endif
+  call premsg('[Cross-section graph]')
+  call premsg('      (section profiles of Mag.Anomaly and Obs.Altitude)')
+  call gparmf('    Specify Width of drawing profiles (in cm) ==> ', wig)
+  call premsg('  [Mag.Anomaly graph scaling]')
+  write(out,'(6x,3(a,f8.1),a)') 'Data Range : [',   &
+    gsmin, ' ,', gsmax, ' ]   Mean :', gssum, ' (nT)'
+  call premsg(out)
+  call gparmf ('    Enter Height of PlotArea (in cm)  ==> ', hig)
+  if (hig <= 0.) call abendm('invalid setting')
+  call gparmf2('    Anom. value range (Bottom and Top) (nT) ==> ', gb, gt)
+  if (gb >= gt) call abendm('invalid setting')
+  call gparma ('      Unit notation for Anom. value ==> ', 8, unit)
+  call premsg('  [Obs.Altitude graph scaling]')
+  write(out,'(6x,3(a,f8.1),a)') 'Data Range : [',   &
+    hsmin, ' ,', hsmax, ' ]   Mean :', hssum, ' (m)'
+  call premsg(out)
+  if (kmd == 1) then
+    ks = 0
+    do i=1,nmd
+      call jcross(mkp(i), pv(1,i), xpa, ypa, xpb, ypb,   &
+                  jj(i), ds1(i), ds2(i), hi1(i), hi2(i))
+!--<input data>
+!   mk:  (model-code)  1 (Block), 2 (hRect), 3 (vLine), 4 (Point)
+!   pv(1-7): [xc, xr, yc, yr, za, zr, vol]  source model locaton
+!   (xpa, ypa) - (xpb, ypb) :
+!     section location of (A)start and (B)end points
+!--<output data>
+!   jj=0: section line is apart from source model
+!   jj=1: section line is close to source model (vLine/Point) at (ds1)
+!          (within small clearance of 1m) in plan view
+!   jj=2: section line intersects or close to source model (Block/hRect)
+!          along line-segment (from ds1 to ds2) in plan view
+!   ds1: distance(1) from (xpa,ypa) (km)
+!   ds2: distance(2) from (xpa,ypa) (km) [0. if jj==1]
+!   hi1: source top altitude (-za) (m)
+!-- hi2: source bottom altitude (-za-zr) (m) [0. if vLine/Point model]
+      if (jj(i) /= 0) cycle
+      deep = hi1(i);  ks = 1;  exit
+    enddo
+    if (ks == 1) then
+      do i=1,nmd
+        call jcross(mkp(i), pv(1,i), xpa, ypa, xpb, ypb,   &
+                    jj(i), ds1(i), ds2(i), hi1(i), hi2(i))
+        if (jj(i) /= 0) then
+          if ((mkp(i)==1) .or. (mkp(i)==3)) then
+            if (hi2(i) < deep) deep = hi2(i)
+          else
+            if (hi1(i) < deep) deep = hi1(i)
+          endif
+        endif
+      enddo
+      write(out,'(6x,a,f8.1,a)')   &
+        'Deepest Source Bottom (along profile) :', deep, ' (m)'
+    endif
+  endif
+  call premsg(out)
+  call premsg('    Enter Height of PlotArea (in cm)')
+  call gparmf('         ( 0 or minus to disappear )  ==> ', hih)
+  if (hih > 0.) then
+    call gparmf2('    Altitude value range (Bottom and Top)  (m) ==> ', hb, ht)
+    if (hb >= ht) call abendm('invalid setting')
+  endif
+  call premsg('')
+  call clspin()
+
+  write(6,'(a)') '----------------------------------------'
+  write(6,'(a)') 'PLXSEC : Plot arbitrary cross-section of anomaly map'
+  write(6,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin))
+  write(6,'(a,a)') ' Output PS filename  : ', fnam(1:lrtrim(fnam))
+  write(6,'(a12,a8,6x,a6,i4,4x,4i8)')   &
+    'Areaname :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
+  write(6,'(14x,2a10,2x,2a8,2a6,a8)')   &
+    'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
+  write(6,'(16x,2f10.2,2i8,2i6,a8)') xs, ys, mszx, mszy, nmx, nmy, salt
+  write(6,'(a)')
+  write(6,'(a,2x,a)') '[Section A-B definition]', sloc
+  write(6,'(4x,a,f9.3,a,f9.3,a,$)') 'A  (', xna, 'N,', yea, 'E )'
+  write(6,'(a,f9.3,a,f9.3,a)') '  -  B  (', xnb, 'N,', yeb, 'E )'
+  write(6,'(2x,a,i5)') 'N0. of points consisting of the section : ', ndp
+  if (kmd == 1) then
+    write(6,'(2x,a)')   &
+      'Source Location is shown, if input data contain such inf.'
+    if (mdl == 0) then
+      write(6,'(4x,a)') '--- CALMAS model inf. not found'
+    else if (mdl == 4) then
+      write(6,'(4x,a)') '--- too many models, truncated to 20'
+    else if (mdl == -1) then
+      write(6,'(4x,a)') '--- invalid CALMAS model inf. neglected'
+    endif
+  endif
+  write(6,'(a)') '[Map-Section Layout / Paper orientation]'
+  write(6,'(4x,4a)') slay, ' layout  on ', shv, ' sheet'
+  write(6,'(a)') '[Index map]'
+  write(6,'(6x,a,i4,4x,a,i4)') 'step=', istp, 'mval=', mval
+  write(6,'(6x,a,f8.0,4x,a3,f8.1,a1,f8.1,a1)')   &
+    'scale  1 /', fscl, shv, him, 'H', wim, 'W'
+  if ((kscl /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
+    write(6,'(6x,a,$)') 'with '
+    if (kshl /= 0) write(6,'(a,$)') '  Coast'
+    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
+      write(6,'(a,f5.1,a,$)') ',', sclu, ' ]'
+    endif
+    write(6,'(a)')
+  endif
+  write(6,'(a)') '[Section graph]'
+  write(6,'(6x,a,f5.1)') 'Width of drawing profiles : ', wig
+  write(6,'(a)') '  [Mag.Anomaly graph scaling]'
+  write(6,'(6x,a,f5.1,$)') 'PlotArea Height:', hig
+  write(6,'(6x,a,f8.1,a,f8.1,a)') 'Bottom/Top value : ', gb, ' /', gt, ' (nT)'
+  if (hih > 0.) then
+    write(6,'(a)') '  [Obs.Altitude graph scaling]'
+    write(6,'(6x,a,f5.1,$)') 'PlotArea Height:', hih
+    write(6,'(6x,a,f8.1,a,f8.1,a)') 'Bottom/Top value : ', hb, ' /', ht, ' (m)'
+  else
+    write(6,'(a)') '  [Obs.Altitude graph is Not shown]'
+  endif
+  write(6,'()')
+  call flush(6)
+!-----------------------------------------------------------------------
+!           (layout [khv])    Horizontal [1]          Vertical [2]
+!     Total height (H:high)   max(him, hig+hih+1)     him+hig+hih+2
+!     Total width  (W:wide)   wim+1+wig               max(wim, wig)
+!  Index map SW corner :      0, (H-him)/2            (W-wim)/2, H-him
+!  Mag.Anom. SW corner :      wim+1, (H-hig+hih+1)/2  (W-wig)/2, hih+1
+!  Obs.Alt. SW corner  :      wim+1, (H-hig-hih-1)/2  (W-wig)/2, 0
+!-----------------------------------------------------------------------
+  if (khv == 1) then
+    wide = wim + 1. + wig;  p1 = 0.;  p2 = wim + 1.;  p3 = p2
+    if (hih > 0.) then
+      high = hih + 1. + hig
+      if (him > high) then
+        q1 = 0.;  q3 = (him-high) / 2.;  q2 = q3 + hih + 1.;  high = him
+      else
+        q1 = (high-him) / 2.;  q3 = 0.;  q2 = q3 + hih + 1.
+      endif
+    else
+      if (him > hig) then
+        high = him;  q1 = 0.;  q2 = (him-hig) / 2.;  q3 = 0.
+      else
+        high = hig;  q1 = (hig-him) / 2.;  q2 = 0.;  q3 = 0.
+      endif
+    endif
+  else
+    if (hih > 0.) then
+      high = hih + 1. + hig + 1. + him
+      wide = wig;  if (wim > wide) wide = wim
+      p1 = (wide-wim) / 2.;  p2 = (wide-wig) / 2.;  p3 = p2
+      q1 = high - him;  q2 = hih + 1.;  q3 = 0
+    else
+      high = hig + 1. + him
+      wide = wig;  if (wim > wide) wide = wim
+      p1 = (wide-wim) / 2.;  p2 = (wide-wig) / 2.;  p3 = p2
+      q1 = high - him;  q2 = 0.;  q3 = 0
+    endif
+  endif
+  call psopn(fnam, shv)
+  call plots(1.5, 1.5)
+  call plot(p1, q1, -3)         !------- Change Origin to [IndexMap]
+  call dfrgbt(0)
+                                        !--- plot IndexMap
+  call dframe(0., 0., wim, him, nmy, nmx)
+  call paintm(g, vlo, vhi, vnul)
+  call newpen(2)
+  call wrect(0., 0., wim, him)
+  call newpen(1)
+  call conts(0., 0., wim, him, nmy, nmx, 0.)
+  call contx(g, vnul, istp, 2, -99999., 99999., 0)
+  if (kscl /= 0) then
+    call plot(sclr,      sclu+0.1, 3)
+    call plot(sclr+dscl, sclu+0.1, 2)
+    call newpen(2)
+    call plot(sclr+dscl, sclu+0.3, 3)
+    call plot(sclr+dscl, sclu,     2)
+    call plot(sclr,      sclu,     2)
+    call plot(sclr,      sclu+0.3, 2)
+    call lstyle('T', 0.5, 0., 0, -255)
+    call ptext(ss, ls, sclr+dscl/2., sclu+0.5, 2)
+    call newpen(1)
+  endif
+  call scisor(0., 0., wim, him)
+  call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
+  if ((klll /= 0) .or. (kshl /= 0)) then
+    xt = xs + float((nmx-1)*mszx)/1000.
+    yt = ys + float((nmy-1)*mszy)/1000.
+    call cviken(float(iorg), float(korg), yo, xo)
+    call cvenik(yo, xs, fi0, fk0)
+    call cvenik(ys, xs, fi1, fk1)
+    call cvenik(yt, xs, fi2, fk2)
+    call cvenik(ys, xt, fi3, fk3)
+    call cvenik(yt, xt, fi4, fk4)
+    call cvenik(yo, xt, fi5, fk5)
+    isf = nint(amin1(fi0,fi1,fi2) - 0.5)
+    itf = nint(amax1(fi3,fi4,fi5) + 0.5)
+    ksf = nint(amin1(fk1,fk3) - 0.5)
+    ktf = nint(amax1(fk2,fk4) + 0.5)
+  endif
+  if (klll /= 0) then
+    do i=isf/iinc,itf/iinc
+      idf = i * iinc
+      if ((idf < isf) .or. (idf > itf)) cycle
+      call cviken(float(idf), float(ksf), y, x)
+      call plot((y-ys)/ccm, (x-xs)/ccm, 3)
+      do kdf=ksf+1,ktf
+        call cviken(float(idf), float(kdf), y, x)
+        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
+      enddo
+    enddo
+    do k=ksf/kinc,ktf/kinc
+      kdf = k * kinc
+      if ((kdf < ksf) .or. (kdf > ktf)) cycle
+      call cviken(float(isf), float(kdf), y, x)
+      call plot((y-ys)/ccm, (x-xs)/ccm, 3)
+      do idf=isf+1,itf
+        call cviken(float(idf), float(kdf), y, x)
+        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
+      enddo
+    enddo
+  endif
+  if (kshl /= 0) then
+    call rshore(mc, isf/60, (itf+59)/60, ksf/60, (ktf+59)/60)
+    call pshore(-ys/ccm, -xs/ccm, kshl, ccm, cviken)
+  endif
+  call scisor(0., 0., 0., 0.)
+                                        !--- plot model location (IndexMap)
+  call newpen(2)
+  do i=1,nmd
+    xc = pv(1,i);  xr = pv(2,i);  yc = pv(3,i);  yr = pv(4,i)
+    za = pv(5,i);  zr = pv(6,i);  k = mkp(i)
+    if ((k==1) .or. (k==2)) then
+      call wrect((yc-yr/2.)/ccm, (xc-xr/2.)/ccm, yr/ccm, xr/ccm)
+    else
+      call wcirc(yc/ccm, xc/ccm, 0.06, 0., 360.)
+    endif
+  enddo
+  call newpen(1)
+!---
+  call penatr(0, 255, 0, 0.06)
+  call plot(yst/ccm, xst/ccm, 3);  call plot(yed/ccm, xed/ccm, 2)
+  call dfpcol(1, 255)
+  call paintc(yst/ccm, xst/ccm, 0.113)
+  call paintr(yed/ccm-0.1, xed/ccm-0.1, 0.2, 0.2)
+
+  xinc = wig / (ndp-1);  pr = ds/wig;  gr = (gt-gb)/hig;  hr = (ht-hb)/hih
+  call plot(p2-p1, q2-q1, -3)   !------- Change Origin to [Mag.Anom profile]
+  call newpen(2)
+                                        !--- plot Mag. profile
+  call wrect(0., 0., wig, hig)
+  call newpen(1)
+  k = 3;  x = 0.
+  do i=1,ndp
+    if (gs(i) == vnul) then
+      k = 3
+    else
+      y = (gs(i)-gb)/gr
+      call plot(x, y, k)
+      k = 2
+    endif
+    x = x + xinc
+  enddo
+  do i=8,1,-1
+    if (unit(i:i) == ' ') cycle
+    n = i + 6; exit
+  enddo
+  call lstyle('T', 0.4, 0., 0, -999)
+  write(buf,'(i6,a8)') nint(gt), unit
+  call ptext(buf, n, wig-1., hig-0.32, 1)
+  write(buf,'(i6,a8)') nint(gb), unit
+  call ptext(buf, n, wig-1., 0.06, 1)
+  call dfpcol(1, 255)
+  call paintc(0., 0., 0.113)
+  call paintr(wig-0.1, -0.1, 0.2, 0.2)
+  call lstyle('T', 0.4, 0., 255, -999)
+  write(buf,'(a1,f0.3,a3,f0.3,a2)') '(', xna, 'N, ', yea, 'E)'
+  call ptext(buf, lrtrim(buf), -0.5, -0.6, 0)
+  write(buf,'(a1,f0.3,a3,f0.3,a2)') '(', xnb, 'N, ', yeb, 'E)'
+  call ptext(buf, lrtrim(buf), wig+0.5, -0.6, 1)
+  write(buf,'(a1,i0,a4)') '[', ndp, 'pts]'
+  call ptext(buf, lrtrim(buf), wig/2., -0.46, 2)
+
+
+  call plot(p3-p2, q3-q2, -3)   !---- Change Origin to [Obs.Alt / Model]
+  if (hih > 0.) then
+    call newpen(2)
+    call wrect(0., 0., wig, hih)
+                                        !--- plot model cross-section
+    if (kmd == 1) then
+      call penatr(0, 255*256, 0, 0.03)
+      do i=1,nmd
+        if (jj(i) == 0) cycle
+        if (jj(i) == 1) then
+          call wcirc(ds1(i)/pr, (hi1(i)-hb)/hr, 0.06, 0., 360.)
+        else
+          call wrect(ds1(i)/pr, (hi2(i)-hb)/hr,   &
+                     (ds2(i)-ds1(i))/pr, (hi1(i)-hi2(i))/hr)
+        endif
+      enddo
+      call newpen(1)
+    endif
+                                        !--- plot Alt. profile
+    k = 3;  x = 0.
+    do i=1,ndp
+      if (hs(i) == hnul) then
+        k = 3
+      else if ((hs(i) < 0.) .and. (alt < 0.)) then
+        k = 3
+      else
+        y = (hs(i)-hb)/hr
+        call plot(x, y, k)
+        k = 2
+      endif
+      x = x + xinc
+    enddo
+    call lstyle('T', 0.4, 0., 0, -999)
+    write(buf,'(i6,a2)') nint(ht), ' m'
+    call ptext(buf, 8, wig-1., hih-0.32, 1)
+    write(buf,'(i6,a2)') nint(hb), ' m'
+    call ptext(buf, 8, wig-1., 0.06, 1)
+    call dfpcol(1, 255)
+    call paintc(0., 0., 0.113)
+    call paintr(wig-0.1, -0.1, 0.2, 0.2)
+  endif
+
+  call plot(-p3, -q3, -3)       !---- Recover Origin to initial state ----
+  call ptext(area, 8, -1., -1., 0)
+  call ptext(fnin, 50, 2., -1., 0)
+  call plote()
+  call pscls()
+
+  if (flog(ldr:ldr) /= ' ') then
+    open(99,file=flog,access='append')
+    write(99,'(//a)') 'PLXSEC : Plot arbitrary cross-section of anomaly map'
+    write(99,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin))
+    write(99,'(a,a)') ' Output PS filename  : ', fnam(1:lrtrim(fnam))
+    write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
+      'Areaname :  ', area, 'Coord.', ncc, iorg,korg,ispa,ispb
+    write(99,'(14x,2a10,2x,2a8,2a6)')   &
+      'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my'
+    write(99,'(16x,2f10.2,2i8,2i6)') xs, ys, mszx, mszy, nmx, nmy
+    write(99,'(a,2x,a)') '[Section A-B definition]', sloc
+    write(99,'(4x,a,f9.3,a,f9.3,a,$)') 'A  (', xna, 'N,', yea, 'E )'
+    write(99,'(a,f9.3,a,f9.3,a)') '  -  B  (', xnb, 'N,', yeb, 'E )'
+    write(99,'(2x,a,i5)') 'No. of points consisting of the section : ', ndp
+    if (kmd == 1) then
+      write(99,'(2x,a)')   &
+        'Source Location is shown, if input data contain such inf.'
+      if (mdl == 0) then
+        write(99,'(4x,a)') '--- CALMAS model inf. not found'
+      else if (mdl == 4) then
+        write(99,'(4x,a)') '--- too many models, truncated to 20'
+      else if (mdl == -1) then
+        write(99,'(4x,a)') '--- invalid CALMAS model inf. neglected'
+      endif
+    endif
+    write(99,'(a)') '[Map-Section Layout / Paper orientation]'
+    write(99,'(4x,4a)') slay, ' layout  on ', shv, ' sheet'
+    write(99,'(a)') '[Index map]'
+    write(99,'(6x,a,i4,4x,a,i4)') 'step=', istp, 'mval=', mval
+    write(99,'(6x,a,f8.0,4x,a3,f8.1,a1,f8.1,a1)')   &
+      'scale  1 /', fscl, shv, him, 'H', wim, 'W'
+    if ((kscl /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
+      write(99,'(6x,a,$)') 'with '
+      if (kshl /= 0) write(99,'(a,$)') '  Coast'
+      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
+        write(99,'(a,f5.1,a,$)') ',', sclu, ' ]'
+      endif
+      write(99,'(a)')
+    endif
+    write(99,'(a)') '[Section graph]'
+    write(99,'(6x,a,f5.1)') 'Width of drawing profiles : ', wig
+    write(99,'(a)') '  [Mag.Anomaly graph scaling]'
+    write(99,'(6x,a,f5.1,$)') 'PlotArea Height:', hig
+    write(99,'(6x,a,f8.1,a,f8.1)') 'Bottom/Top value : ', gb, ' /', gt
+    if (hih > 0.) then
+      write(99,'(a)') '  [Obs.Altitude graph scaling]'
+      write(99,'(6x,a,f5.1,$)') 'PlotArea Height:', hih
+      write(99,'(6x,a,f8.1,a,f8.1)') 'Bottom/Top value : ', hb, ' /', ht
+    else
+      write(99,'(a)') '  [Obs.Altitude graph is Not shown]'
+    endif
+    call strdtm(sdtm)
+    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
+    close(99)
+  endif
+  stop
+  end
+!=======================================================================
+!---
+  subroutine jcross(mk, pv, xpa, ypa, xpb, ypb, jj, ds1, ds2, hi1, hi2)
+  implicit none
+  integer,intent(in) :: mk
+  real,intent(in) :: pv(7)
+  real,intent(in) :: xpa, ypa, xpb, ypb
+  integer,intent(out) :: jj
+  real,intent(out) :: ds1, ds2, hi1, hi2
+!--<input data>
+!   mk:  (model-code)  1 (Block), 2 (hRect), 3 (vLine), 4 (Point)
+!   pv(1-7): [xc, xr, yc, yr, za, zr, vol]  source model locaton
+!   (xpa, ypa) - (xpb, ypb) :
+!     section location of (A)start and (B)end points
+!--<output data>
+!   jj=0: section line is apart from source model
+!   jj=1: section line is close to source model (vLine/Point) at (ds1)
+!          (within small clearance of 1m) in plan view
+!   jj=2: section line intersects or close to source model (Block/hRect)
+!          along line-segment (from ds1 to ds2) in plan view
+!   ds1: distance(1) from (xpa,ypa) (km)
+!   ds2: distance(2) from (xpa,ypa) (km) [0. if jj==1]
+!   hi1: source top altitude (-za) (m)
+!-- hi2: source bottom altitude (-za-zr) (m) [0. if vLine/Point model]
+  integer :: kx, j
+  real :: dx, dy, ds, rd1, rd2, xp1, yp1, xp2, yp2
+  real :: sxl, sxh, syl, syh, x0, y0, dist, d0pa, d0pb
+  dx = xpb-xpa;  dy = ypb-ypa;  ds = sqrt(dx*dx + dy*dy)
+  kx = 0
+  if ((mk==3) .or. (mk==4)) then        !-- vLine/Point model
+!---------------------------------------------------------------------------
+!   General formula of the Line passing two points (xpa,ypa) and (xpb,ypb)
+!    is expressed as  [ Ax + By + C = 0 ] ,
+!    where  A = ypa - ypb , B = xpb - xpa , C = xpa*ypb - xpb*ypa .
+!   The distance from a point (x0,y0) to this line is 
+!      |A x0�{B y0�{C| / sqrt(A*A + B*B)
+!   A ==> (-dy), B ==> dx, sqrt(A*A + B*B) ==> ds �ƑΉ�����D
+!---------------------------------------------------------------------------
+    x0 = pv(1);  y0 = pv(3)
+    dist = abs((xpa*ypb-xpb*ypa) - dy*x0 + dx*y0) / ds
+    if (dist < 0.001) then
+      d0pa = sqrt((x0-xpa)*(x0-xpa) + (y0-ypa)*(y0-ypa) - dist*dist)
+      d0pb = sqrt((x0-xpb)*(x0-xpb) + (y0-ypb)*(y0-ypb) - dist*dist)
+      if (d0pa < d0pb)then
+        if (d0pb <= ds) then
+          kx = 1;  ds1 = d0pa;
+        endif
+      else
+        if (d0pa <= ds) then
+          kx = 1;  ds1 = ds - d0pb;
+        endif
+      endif
+    endif
+    if (kx == 1) then
+      jj = 1;  ds2 = 0.;  hi1 = -pv(5)*1000.;  hi2 = 0.
+    else 
+      jj = 0
+    endif
+  else                                  !-- Block/hRect model
+    sxl = pv(1) - pv(2)/2.;  sxh = pv(1) + pv(2)/2.
+    syl = pv(3) - pv(4)/2.;  syh = pv(3) + pv(4)/2.
+    if (abs(dx) > abs(dy)) then         !---- section nearer to N-S
+      rd1 = (sxl-xpa)/dx;  rd2 = (sxh-xpa)/dx
+      if ((rd1>1.) .and. (rd2>1.)) then
+      else if ((rd1<0.) .and. (rd2<0.)) then
+      else
+        yp1 = ypa + dy*rd1;  yp2 = ypa + dy*rd2
+        if ((yp1>syh) .and. (yp2>syh)) then
+        else if ((yp1<syl) .and. (yp2<syl)) then
+        else
+          if (yp1 > syh) then
+            rd1 = (syh-ypa)/dy; xp1 = xpa + dx*rd1; yp1 = ypa + dy*rd1
+          endif
+          if (yp2 > syh) then
+            rd2 = (syh-ypa)/dy; xp2 = xpa + dx*rd2; yp2 = ypa + dy*rd2
+          endif
+          if (yp1 < syl) then
+            rd1 = (syl-ypa)/dy; xp1 = xpa + dx*rd1; yp1 = ypa + dy*rd1
+          endif
+          if (yp2 < syl) then
+            rd2 = (syl-ypa)/dy; xp2 = xpa + dx*rd2; yp2 = ypa + dy*rd2
+          endif
+          kx = 2
+        endif
+      endif
+    else                                !---- section nearer to E-W
+      rd1 = (syl-ypa)/dy;  rd2 = (syh-ypa)/dy
+      if ((rd1>1.) .and. (rd2>1.)) then
+      else if ((rd1<0.) .and. (rd2<0.)) then
+      else
+        xp1 = xpa + dx*rd1;  xp2 = xpa + dx*rd2
+        if ((xp1>sxh) .and. (xp2>sxh)) then
+        else if ((xp1<sxl) .and. (xp2<sxl)) then
+        else
+          if (xp1 > sxh) then
+            rd1 = (sxh-xpa)/dx; yp1 = ypa + dy*rd1; xp1 = xpa + dx*rd1
+          endif
+          if (xp2 > sxh) then
+            rd2 = (sxh-xpa)/dx; yp2 = ypa + dy*rd2; xp2 = xpa + dx*rd2
+          endif
+          if (xp1 < sxl) then
+            rd1 = (sxl-xpa)/dx; yp1 = ypa + dy*rd1; xp1 = xpa + dx*rd1
+          endif
+          if (xp2 < sxl) then
+            rd2 = (sxl-xpa)/dx; yp2 = ypa + dy*rd2; xp2 = xpa + dx*rd2
+          endif
+          kx = 2
+        endif
+      endif
+    endif
+    if (kx == 2) then
+      jj = 2;  hi1 = -pv(5)*1000.;  hi2 = -(pv(5)+pv(6))*1000.
+      if (rd1 < rd2) then
+        ds1 = ds*rd1; ds2 = ds*rd2
+      else
+        ds1 = ds*rd2; ds2 = ds*rd1
+      endif
+    else
+      do j=1,4
+        if (j < 3) then; x0 = sxl; else; x0 = sxh; endif
+        if ((j-j/2*2) == 1) then; y0 = syl; else; y0 = syh; endif
+        dist = abs((xpa*ypb-xpb*ypa) - dy*x0 + dx*y0) / ds
+        if (dist < 0.001) then
+          d0pa = sqrt((x0-xpa)*(x0-xpa) + (y0-ypa)*(y0-ypa) - dist*dist)
+          d0pb = sqrt((x0-xpb)*(x0-xpb) + (y0-ypb)*(y0-ypb) - dist*dist)
+          if (d0pa < d0pb) then
+            if (d0pb <= ds) then
+              if (kx == 0) then
+                kx = 1; ds1 = d0pa; ds2 = 0.
+              else
+                kx = 2; ds2 = d0pa; exit
+              endif
+            endif
+          else
+            if (d0pa <= ds) then
+              if (kx == 0) then
+                kx = 1; ds1 = ds - d0pb; ds2 = 0.
+              else
+                kx = 2; ds2 = ds - d0pb; exit
+              endif
+            endif
+          endif
+        endif
+      enddo
+      if (kx > 0) then
+        jj = kx;  hi1 = -pv(5)*1000.;  hi2 = 0.
+      else
+        jj = 0
+      endif
+    endif
+  endif
+  return
+  end subroutine
diff -urN ORG/gdmp/shade.f90 gdmp/shade.f90
--- ORG/gdmp/shade.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/shade.f90	2019-12-17 10:44:00.000000000 +0900
@@ -5,9 +5,9 @@
   real,parameter :: RAD = (180./3.141593)
   character(80) :: wdr, flog, fnin, fnam, buf
   character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
-  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
+  character(6) :: ss;  character(3) :: spl = '   ';  character(1) :: yn
   integer :: ldr, i, j, n, k, IER
-  integer :: kpl, kscl, klll, kshl, km, iinc, kinc
+  integer :: kpl, kscl, klll, kshl, km, ls, iinc, kinc
   integer :: nc, ncc, mc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, mx, my
   integer :: isf, itf, ksf, ktf, idf, kdf
@@ -132,17 +132,55 @@
       fscl = fscl * 100.
       ccm = fscl/100000.
       kscl = 0;  sclr = 0.;  sclu = 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
@@ -257,9 +295,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)
diff -urN ORG/gdmp/xplmap.f90 gdmp/xplmap.f90
--- ORG/gdmp/xplmap.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/xplmap.f90	2019-12-17 10:44:00.000000000 +0900
@@ -10,12 +10,12 @@
   character(LBUF) :: sbf
   character(80) :: wdr, flog, fnin, fnobj, fncap, fnam, buf
   character(72) :: out;  character(40) :: slbl;  character(20) :: sdtm
-  character(8) :: area;  character(4) :: sht, s4;  character(3) :: fnt
-  character(1) :: yn
+  character(8) :: area;  character(6) :: ss;  character(4) :: sht
+  character(3) :: fnt;  character(1) :: yn
   integer :: ldr, nseq, ns, n, i, j, k, kpolyg, m, nct, IER
   integer :: nc, ncc, mc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
-  integer :: istep, istp, lcapt, km, iinc, kinc, kscl, klll, kshl, klbl
+  integer :: istep, istp, lcapt, km, ls, iinc, kinc, kscl, klll, kshl, klbl
   integer :: ibcol, itcol, idf, isf, itf, kdf, ksf, ktf, ksp, mark
   real :: vnul, dmy, xs, ys, gsum, gmin, gmax, szx, szy, swd, shi
   real :: fscl, am, bm, csiz, high, wide, rs, ts, ccm, dscl, sclr, sclu
@@ -171,21 +171,55 @@
   call premsg(out)
 
   ccm = fscl/100000.
-  km = nint(fscl)/25000
-  if (km < 1) km = 1
-  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
-  if ((km >= 50) .and. (km < 120)) km = 60
-  if ((km > 120) .and. (km < 300)) km = 120
-  if ((km > 300) .and. (km < 600)) km = 300
-  if (km > 600) km = 600
-  iinc = km
-  kinc = km + km/2
-  if (km > 50) km = 50
-  if ((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
 
   kscl = 0;  sclr = 0.;  sclu = 0.;  klll = 0;  kshl = 0
   call prompt('  Draw  Scale Bar  ?  ')
@@ -248,8 +282,10 @@
     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) write(6,'(a,f5.1,a,f5.1,a,$)')   &
-                     '  ScaleBar [', sclr, ',', sclu, ' ]'
+    if (kscl /= 0) then
+      write(6,'(a,f5.1,$)') '  ScaleBar [', sclr
+      write(6,'(a,f5.1,a,$)') ',', sclu, ' ]'
+    endif
     write(6,'(a)')
   endif
   if (klbl /= 0) write(6,'(a,a)') '      Identifier : ', slbl
@@ -415,9 +451,8 @@
 
   if (kscl /= 0) then
     call premsg('--- Drawing Scale-Bar')
-    write(s4,'(i2,a2)') km, 'km'
     call lstyle('TB ', 0.7, 0., 0, -999)
-    call ptext(s4, 4, sclr+dscl/2., sclu+0.5, 2)
+    call ptext(ss, ls, sclr+dscl/2., sclu+0.5, 2)
     call newpen(2)
     call plot(sclr,      sclu+0.1, 3)
     call plot(sclr+dscl, sclu+0.1, 2)
@@ -579,8 +614,10 @@
       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) write(99,'(a,f5.1,a,f5.1,a,$)')   &
-                       '  ScaleBar [', sclr, ',', sclu, ' ]'
+      if (kscl /= 0) then
+        write(99,'(a,f5.1,$)') '  ScaleBar [', sclr
+        write(99,'(a,f5.1,a,$)') ',', sclu, ' ]'
+      endif
       write(99,'(a)')
     endif
     if (klbl /= 0) write(99,'(a,a)') '      Identifier : ', slbl
diff -urN ORG/gdmp/xplmapc.f90 gdmp/xplmapc.f90
--- ORG/gdmp/xplmapc.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/xplmapc.f90	2019-12-17 10:44:00.000000000 +0900
@@ -10,12 +10,12 @@
   character(LBUF) :: sbf
   character(80) :: wdr, flog, fnin, fnobj, fncap, fnam, buf
   character(72) :: out;  character(40) :: slbl;  character(20) :: sdtm
-  character(10) :: cunit;  character(8) :: area;  character(6) :: s6
+  character(10) :: cunit;  character(8) :: area;  character(6) :: ss
   character(4) :: sht;  character(3) :: fnt;  character(1) :: yn
   integer :: ldr, nseq, ns, n, i, j, k, kpolyg, m, nct, IER
   integer :: nc, ncc, mc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
-  integer :: istep, istp, km, iinc, kinc, kscl, klll, kshl, klbl
+  integer :: istep, istp, km, ls, iinc, kinc, kscl, klll, kshl, klbl
   integer :: ibcol, itcol, idf, isf, itf, kdf, ksf, ktf, ksp, mark
   integer :: kcgl, mval, nras
   real :: vnul, dmy, xs, ys, gsum, gmin, gmax, szx, szy, swd, shi
@@ -178,21 +178,55 @@
   call premsg(out)
 
   ccm = fscl/100000.
-  km = nint(fscl)/25000
-  if (km < 1) km = 1
-  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
-  if ((km >= 50) .and. (km < 120)) km = 60
-  if ((km > 120) .and. (km < 300)) km = 120
-  if ((km > 300) .and. (km < 600)) km = 300
-  if (km > 600) km = 600
-  iinc = km
-  kinc = km + km/2
-  if (km > 50) km = 50
-  if ((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
 
   kcgl = 0;  cglf = 1.;  cglr = 0.;  cglu = 0.
   call prompt('  Draw Color-Grading Legend ?')
@@ -274,8 +308,10 @@
     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) write(6,'(a,f5.1,a,f5.1,a,$)')   &
-                     '  ScaleBar [', sclr, ',', sclu, ' ]'
+    if (kscl /= 0) then
+      write(6,'(a,f5.1,$)') '  ScaleBar [', sclr
+      write(6,'(a,f5.1,a,$)') ',', sclu, ' ]'
+    endif
     write(6,'(a)')
   endif
   if (klbl /= 0) write(6,'(a,a)') '      Identifier : ', slbl
@@ -445,9 +481,8 @@
 
   if (kscl /= 0) then
     call premsg('--- Drawing Scale-Bar')
-    write(s6,'(i2,a2)') km, 'km'
     call lstyle('TB ', 0.7, 0., 0, -999)
-    call ptext(s6, 4, sclr+dscl/2., sclu+0.5, 2)
+    call ptext(ss, ls, sclr+dscl/2., sclu+0.5, 2)
     call newpen(2)
     call plot(sclr,      sclu+0.1, 3)
     call plot(sclr+dscl, sclu+0.1, 2)
@@ -524,8 +559,8 @@
     call lstyle('CFN', cglf*0.48, 0., 0, -999)
     do i=-10,10
       y = y + cglf
-      write(s6,'(i6)') mval+istp*2*i
-      call ptext(s6, 6, x-cglf*0.06, y-cglf*0.12, 1)
+      write(ss,'(i6)') mval+istp*2*i
+      call ptext(ss, 6, x-cglf*0.06, y-cglf*0.12, 1)
     enddo
   endif
 
@@ -641,8 +676,10 @@
       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) write(99,'(a,f5.1,a,f5.1,a,$)')   &
-                       '  ScaleBar [', sclr, ',', sclu, ' ]'
+      if (kscl /= 0) then
+        write(99,'(a,f5.1,$)') '  ScaleBar [', sclr
+        write(99,'(a,f5.1,a,$)') ',', sclu, ' ]'
+      endif
       write(99,'(a)')
     endif
     if (klbl /= 0) write(99,'(a,a)') '      Identifier : ', slbl
diff -urN ORG/gdmp/xplmapcs.f90 gdmp/xplmapcs.f90
--- ORG/gdmp/xplmapcs.f90	2017-12-23 18:00:00.000000000 +0900
+++ gdmp/xplmapcs.f90	2019-12-17 10:44:00.000000000 +0900
@@ -10,12 +10,12 @@
   character(LBUF) :: sbf
   character(80) :: wdr, flog, fnin, fnobj, fncap, fnam, buf
   character(72) :: out;  character(40) :: slbl;  character(20) :: sdtm
-  character(10) :: cunit;  character(8) :: area;  character(6) :: s6
+  character(10) :: cunit;  character(8) :: area;  character(6) :: ss
   character(4) :: sht;  character(3) :: fnt;  character(1) :: yn
   integer :: ldr, nseq, ns, n, i, j, k, kpolyg, m, nct, IER
   integer :: nc, ncc, mc, iorg, korg, ispa, ispb
   integer :: ixs, iys, mszx, mszy, nmx, nmy
-  integer :: istep, istp, km, iinc, kinc, kscl, klll, kshl, klbl
+  integer :: istep, istp, km, ls, iinc, kinc, kscl, klll, kshl, klbl
   integer :: ibcol, itcol, idf, isf, itf, kdf, ksf, ktf, ksp, mark
   integer :: kcgl, mval, nras
   real :: vnul, dmy, xs, ys, gsum, gmin, gmax, szx, szy, swd, shi
@@ -189,21 +189,55 @@
   call premsg(out)
 
   ccm = fscl/100000.
-  km = nint(fscl)/25000
-  if (km < 1) km = 1
-  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
-  if ((km >= 50) .and. (km < 120)) km = 60
-  if ((km > 120) .and. (km < 300)) km = 120
-  if ((km > 300) .and. (km < 600)) km = 300
-  if (km > 600) km = 600
-  iinc = km
-  kinc = km + km/2
-  if (km > 50) km = 50
-  if ((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
 
   kcgl = 0;  cglf = 1.;  cglr = 0.;  cglu = 0.
   call prompt('  Draw Color-Shading Legend ?')
@@ -288,8 +322,10 @@
     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) write(6,'(a,f5.1,a,f5.1,a,$)')   &
-                     '  ScaleBar [', sclr, ',', sclu, ' ]'
+    if (kscl /= 0) then
+      write(6,'(a,f5.1,$)') '  ScaleBar [', sclr
+      write(6,'(a,f5.1,a,$)') ',', sclu, ' ]'
+    endif
     write(6,'(a)')
   endif
   if (klbl /= 0) write(6,'(a,a)') '      Identifier : ', slbl
@@ -500,9 +536,8 @@
 
   if (kscl /= 0) then
     call premsg('--- Drawing Scale-Bar')
-    write(s6,'(i2,a2)') km, 'km'
     call lstyle('TB ', 0.7, 0., 0, -999)
-    call ptext(s6, 4, sclr+dscl/2., sclu+0.5, 2)
+    call ptext(ss, ls, sclr+dscl/2., sclu+0.5, 2)
     call newpen(2)
     call plot(sclr,      sclu+0.1, 3)
     call plot(sclr+dscl, sclu+0.1, 2)
@@ -578,8 +613,8 @@
     call lstyle('CFN', cglf*0.48, 0., 0, -999)
     do i=-10,10
       y = y + cglf
-      write(s6,'(i6)') mval+istp*2*i
-      call ptext(s6, 6, x-cglf*0.06, y-cglf*0.12, 1)
+      write(ss,'(i6)') mval+istp*2*i
+      call ptext(ss, 6, x-cglf*0.06, y-cglf*0.12, 1)
     enddo
   endif
 
@@ -698,8 +733,10 @@
       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) write(99,'(a,f5.1,a,f5.1,a,$)')   &
-                       '  ScaleBar [', sclr, ',', sclu, ' ]'
+      if (kscl /= 0) then
+        write(99,'(a,f5.1,$)') '  ScaleBar [', sclr
+        write(99,'(a,f5.1,a,$)') ',', sclu, ' ]'
+      endif
       write(99,'(a)')
     endif
     if (klbl /= 0) write(99,'(a,a)') '      Identifier : ', slbl
