program ch3303
use dislin
logical :: trial, screen
real :: long, lat

! switch diagnostics  on or off. this was used when developing
! the original uniras version of the code. i haven't needed
! to switch it on with the dislin version.

  screen = .false.
  trial = .false.

!      print *,' testing on/off'
!      read (unit=*, fmt=*, end=9,err=9)trial
!9     continue

!     read in the tsunami data

  call datain(trial)

!     i now have all the tsunami data latitude and longitude
!     values read in to the arrays in the tsunam common block.


! the following options were available with the original
! uniras version of the code. they were not available with
! the dislin version. they are therefore commented out.
!     the next thing to do is sort out what kind of
!     map projection the user wants. i am offering
!     one of the following
!     projection
!     lambert       - equal area        -  rectangle
!     mercator      - equal direction   -  rectangle
!     hammer        - equal area        -  oval
!     bonne         -                   -  heart
!     orthographic  - globe             -  round
!  print *, ' what projection would you like?'
!  print *
!  print *, ' 1 = lambert       - equal area        -  rectangle'
!  print *, ' 2 = mercator      - equal direction   -  rectangle'
!  print *, ' 3 = hammer        - equal area        -  oval'
!  print *, ' 4 = bonne         -                   -  heart'
!  print *, ' 5 = orthographic  - globe             -  round'
!100 read (unit=*,fmt=*,end=110,err=110) iproj
!110 if ((iproj<1) .or. (iproj>5)) then
!    print *, ' please input a number in the range 1 to 5'
!    go to 100
!  end if
!     the next thing to do is set the centre of the map
!     in this case the pacific

  iproj=1
  lat = 0.0
  long = 180.0

! the following options were available with the original
! uniras version of the code. they were not available with
! the dislin version. they are therefore commented out.
!     i offer the user a choice of region
!     to plot.
!  print *, ' which region do you wish to plot?'
!
!  print *, ' 1 = hawaii'
!  print *, ' 2 = new zealand and south pacific islands'
!  print *, ' 3 = papua new guinea and solomon islands'
!  print *, ' 4 = indonesia'
!  print *, ' 5 = philippines'
!  print *, ' 6 = japan'
!  print *, ' 7 = kuril islands and kamchatka'
!
!  print *, ' 9 = west coast - north and central america'
!  print *, ' 10 = west coast - south america'
!120 read (unit=*,fmt=*,end=130,err=130) nreg
!130 if ((nreg<0) .or. (nreg>10)) then
!    print *, ' please input a number between 0 and 10 inclusive'
!    go to 120
!  end if

  nreg=0

! dislin initialisation routines and setting of some basic components
! of the plot. these are based on two sample dislin programs.

! initialise dislin

  call disini

! choose font

  call psfont('times-roman')

! determines the position of an axis system.
! the lower left corner of the axis system

  call axspos(400,1850)

! the size of the axis system
! are the length and height of an axis system in plot coordinates. the
default
! values are set to 2/3 of the page length and height.

  call axslen(2400,1400)

! define axis title

  call name('longitude','x')

! define axis title

  call name('latitude','y')

! this routine plots a title over an axis system.

  call titlin('plot of 3034 tsunami events ',3)

! determines which label types will be plotted on an axis.
! map defines geographical labels which are plotted as non negative
floating-point
! numbers with the following characters 'w', 'e', 'n' and 's'.

  call labels('map','xy')

! plots a geographical axis system.

  call grafmp(-180.,180.,-180.,90.,-90.,90.,-90.,30.)

! the statement call gridmp (i, j) overlays an axis system with a longitude
! and latitude grid where i and j are the number of grid lines between
labels in
! the x- and y-direction.

  call gridmp(1,1)

! the routine world plots coastlines and lakes.

  call world

! the angle and height of the characters can be changed with the routines
! angle and height.

  call height(50)

! this routine plots a title over an axis system.
! the title may contain up to four lines of text designated
! with titlin.

  call title

! this is a call to the convert routine.
! this was required by uniras
!  call convrt(trial)

! this is a call to the routine that actually plots each event.

  call plotem(trial,nreg)

! disfin terminates dislin and prints a message on the screen. the level is
set back to 0.

  call disfin

end program ch3302

subroutine datain(trial)
 common /tsunam/ &
  reg0la(378) , &
  reg0lo(378) , &
  reg1la(206) , &
  reg1lo(206) , &
  reg2la(41) , &
  reg2lo(41) , &
  reg3la(54) , &
  reg3lo(54) , &
  reg4la(60) , &
  reg4lo(60) , &
  reg5la(1540) , &
  reg5lo(1540) , &
  reg6la(80) , &
  reg6lo(80) , &
  reg7la(144) , &
  reg7lo(144) , &
  reg8la(245) , &
  reg8lo(245) , &
  reg9la(285) , &
  reg9lo(285)

! this subroutine reads in the tsunami data
! uniras uses the following unit numbers
! 5,6,7,20,21,22,24,25,26,27,28,33
! so i have used 50.

logical :: trial
character (80) :: filnam
  if (trial) then
    print *, ' entering data input phase'
  end if
  filnam = 'tsunami.dat'
  open (unit=50,file=filnam,err=100,status='old')
  go to 110
100 print *, ' error opening data file'
  print *, ' program terminates'
  stop
110 do i = 1, 378
    read (unit=50,fmt=1000) reg0la(i), reg0lo(i)
  end do
1000 format (1x,f7.2,2x,f7.2)
  do i = 1, 206
    read (unit=50,fmt=1000) reg1la(i), reg1lo(i)
  end do
  do i = 1, 41
    read (unit=50,fmt=1000) reg2la(i), reg2lo(i)
  end do
  do i = 1, 54
    read (unit=50,fmt=1000) reg3la(i), reg3lo(i)
  end do
  do i = 1, 60
    read (unit=50,fmt=1000) reg4la(i), reg4lo(i)
  end do
  do i = 1, 1540
    read (unit=50,fmt=1000) reg5la(i), reg5lo(i)
  end do
  do i = 1, 80
    read (unit=50,fmt=1000) reg6la(i), reg6lo(i)
  end do
  do i = 1, 144
    read (unit=50,fmt=1000) reg7la(i), reg7lo(i)
  end do
  do i = 1, 245
    read (unit=50,fmt=1000) reg8la(i), reg8lo(i)
  end do
  do i = 1, 285
    read (unit=50,fmt=1000) reg9la(i), reg9lo(i)
  end do
  if (trial) then
    do i = 1, 10
      print *, reg0la(i), '   ', reg0lo(i)
    end do
    print *, ' exiting data input phase'
    read *, dummy
  end if
end subroutine datain

subroutine plotem(trial,nreg)
use dislin
 common /tsunam/ &
  reg0la(378) , &
  reg0lo(378) , &
  reg1la(206) , &
  reg1lo(206) , &
  reg2la(41) , &
  reg2lo(41) , &
  reg3la(54) , &
  reg3lo(54) , &
  reg4la(60) , &
  reg4lo(60) , &
  reg5la(1540) , &
  reg5lo(1540) , &
  reg6la(80) , &
  reg6lo(80) , &
  reg7la(144) , &
  reg7lo(144) , &
  reg8la(245) , &
  reg8lo(245) , &
  reg9la(285) , &
  reg9lo(285)

! this subroutine plots all of the tsunamis onto the map as coloured
! points, with a different colour per region. i have chosen
! a dot size of 1 mm, and step through the colour pallette.
! the default may not be appropriate.

  logical :: trial
  integer :: nreg
  integer :: kolour=10
  data dwidth/1.0/
  if (trial) then
    dwidth = 5.0
    print *, ' entering plot points'
  end if
  call incmrk(-1)
  if (nreg==0) then
    call setclr(kolour)
    call curvmp(reg0lo,reg0la,378)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg1lo,reg1la,206)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg2lo,reg2la,41)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg3lo,reg3la,54)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg4lo,reg4la,60)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg5lo,reg5la,1540)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg6lo,reg6la,80)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg7lo,reg7la,144)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg8lo,reg8la,245)
    kolour = kolour +30
    call setclr(kolour)
    call curvmp(reg9lo,reg9la,285)
  else if (nreg==1) then
    kolour = 10
    call setclr(kolour)
    call curvmp(reg0lo,reg0la,378)
  else if (nreg==2) then
    kolour = 20
    call setclr(kolour)
    call curvmp(reg1lo,reg1la,206)
  else if (nreg==3) then
    kolour = 30
    call setclr(kolour)
    call curvmp(reg2lo,reg2la,41)
  else if (nreg==4) then
    kolour = 40
    call setclr(kolour)
    call curvmp(reg3lo,reg3la,54)
  else if (nreg==5) then
    kolour = 50
    call setclr(kolour)
    call curvmp(reg4lo,reg4la,60)
  else if (nreg==6) then
    kolour = 60
    call setclr(kolour)
    call curvmp(reg5lo,reg5la,1540)
  else if (nreg==7) then
    kolour = 70
    call setclr(kolour)
    call curvmp(reg6lo,reg6la,80)
  else if (nreg==8) then
    kolour = 80
    call setclr(kolour)
    call curvmp(reg7lo,reg7la,144)
  else if (nreg==9) then
    kolour = 90
    call setclr(kolour)
    call curvmp(reg8lo,reg8la,245)
  else if (nreg==10) then
    kolour = 100
    call setclr(kolour)
    call curvmp(reg9lo,reg9la,285)
  end if
  if (trial) then
    print *, ' exiting plot points'
  end if

end subroutine plotem