! This version verified 12th Nov 2003 
!                                                * pfig.f90    *
!                                                * Purser 2001 *
!                    Last major revision:  28th July     2000
!                    Last minor revisions: 20th December 2002
!                                          24th June     2003

module pfkinds
  integer, parameter :: sp = kind(1.0)
  integer, parameter :: dp = kind(1.0d0)
end module pfkinds

!==============================================================================

module pfconsts  
  use pfkinds
  implicit none
  integer,parameter:: npen=20,ntnr=20,nltp=20,nmtp=20,&
       nmtype=10,nltype=20,nfonts=7,mtypei_def=0,iframe_max=99,jv_max=200,&
       mpeni_def=0,ltypei_def=0,lpeni_def=0,ggunit_def=1,&
       ncon_max=100,npatch=32,npatchp=npatch+1

  real(sp),parameter:: dy=8128., c30=.8660254, s30=.5,&
       scaleh=99.096576,scalew=182.88,scalewh=91.44,chh_def=1.,chw_def=.5,&
       xrun_def=1.,yrise_def=0.,&
       xc2_def=0.,yc2_def=0.,&                       ! Pen position
       xw0_def=0.,xw1_def=1.,yw0_def=0.,yw1_def=1.,& ! Window bounds
       xv0_def=0.,xv1_def=1.,yv0_def=0.,yv1_def=1.,& ! Viewport bounds
       wid_def=.5,cfudge=1.2192

  real(dp),parameter:: c30_d=.866025403784439d0,s30_d=.5d0

  character*4,parameter:: fistem_def='figr'
  logical,parameter:: lgmess_def=.true.
end module pfconsts

module pfstatus !=============================================================
  use pfkinds
  use pfconsts
  implicit none
  integer:: jv,iframe,ggunit,&
       ifonts,ifonts_a,&                    ![1,7]
       kfine,kfine_a, &                     ![1,5]
       imtype,iltype,iltype_old,&           ![0,nmtype],[0,nltype],[0,nltype]
       impen,ilpen,itpen,igpen,ispen,ilpen_old,itpen_old,&    ![0,npen]
       itnr,iltp,imtp,itick,&               ![0,ntnr],[0,nltp],[0,nmtp],[-1,1]
       nbdiv,nbdivp,&
       modec,mcon,ncon,&                    ![0,1],[1,ncon_max],[1,ncon_max]
       ipos1,ineg1,ipos5,ineg5,&            ![0,nltp]
       ksty,&                               ![1,2]
       tmx,tmy,tnx,tny,tnxm,tnym,mnnx
  integer,dimension(ncon_max)       :: ctyp         ![0,nltp]
  integer,dimension(npatchp,npatchp):: kt
  integer,dimension(npatch)         :: kty
  integer,dimension(0:nmtp)         :: mtypei,mpeni ![0,nmtype],[0,npen]
  integer,dimension(0:nltp)         :: ltypei,lpeni ![0,nltype],[0,npen]
  integer,dimension(0:nltype)       :: ndtlen       ![1,3]
  real(sp):: ax,ay,bx,by,xc2,yc2,xc2_a,yc2_a,&
       chh,chw,dchh,dchw,chh_a,xrun,yrise,rpli,rpmi,&
       wid_a,hue_a,sat_a,bri_a,&
       xcapb,ycapb, xcapl,ycapl, xcapt,ycapt, xcapr,ycapr,&
       xw0_old,xw1_old,yw0_old,yw1_old,&
       x0cwn,x1cwn,y0cwn,y1cwn,&
       rsdar,rsodd,&
       dinck,dincj,acon1,acon4,sx,sy,px,py
  real(sp),dimension(2:4,4)     :: fw1
  real(sp),dimension(3,4)       :: fwn
  real(sp),dimension(4,4)       :: fw     
  real(sp),dimension(0:ntnr)    :: xw0,xw1,yw0,yw1,  xv0,xv1,yv0,yv1
  real(sp),dimension(0:npen)    :: pwid,phue,psat,pbri
  real(sp),dimension(3,0:nltype):: dtlen1,dtlen2
  real(sp),dimension(ncon_max)  :: cval
  real(sp),dimension(jv_max)    :: xarr,yarr,zarr
  character*4 fistem
  character*9 finame
  logical lini,lopen,lgmess,lcolor,lstroke,llabel,lcwn,ltform,lcval,lportrait
  data lini/.false./,lportrait/.false./

end module pfstatus !=========================================================
 
module pfig
  private ggjplot,ggywtod,ggxwtod,ggxdtow,ggydtow,ggyttow,ggxttow, &
       ggscaleh,ggscalew,ggdchh,ggset_status,ggscapb,ggscapt,ggscapl,ggscapr,&
       ggoappen,ggappen,ggcurve,ggkrcon,gglrcon,ggindef_a,ggopen,ggsquash,&
       ggdecade,ggpl2,gghpl,ggdefaults,ggenlin,ggoenlin,ggupdate,ggupdatex,&
       gginixax,ggfinxax,gglline2,gglline3,ggfline1,ggfline2,ggfline3,&
       ggsetw,ggshade,ggscapxy,gguform,ggfmt99,ggfmt100,ggfmt101,ggfmt102,&
       ggfmt103,ggfmt104,ggfmt105,ggfmt106,ggfmt107,ggfmt108,ggfmt109,&
       ggfmt112,ggfmt127,ggfmt199,ggjwrite,ggjclose,&
       hpl2,htpl2,hpl3,gtcon,gopm,grect0,grect4,grange1,grange2,gtrange,&
       gtxax,gtyax
  interface ggxwtod
     module procedure ggxwtod
  end interface
  interface ggywtod
     module procedure ggywtod
  end interface
  interface ggxdtow
     module procedure ggxdtow
  end interface
  interface ggydtow
     module procedure ggydtow
  end interface
  interface ggyttow
     module procedure ggyttow
  end interface
  interface ggxttow
     module procedure ggxttow
  end interface
  interface ggscaleh
     module procedure ggscaleh
  end interface
  interface ggscalew
     module procedure ggscalew
  end interface
  interface ggdchh
     module procedure ggdchh
  end interface
  interface wstring
     module procedure wstring
  end interface
  interface gpink
     module procedure gpink
  end interface
  interface gblue
     module procedure gblue
  end interface
  interface gorange
     module procedure gorange
  end interface
  interface ggreen
     module procedure ggreen
  end interface
  interface gpurple
     module procedure gpurple
  end interface
  interface gred
     module procedure gred
  end interface

  interface gpl
     module procedure hpl2,htpl2,hpl3,hpl2_d,htpl2_d,hpl3_d
  end interface
  interface gpm
     module procedure gpm,gopm,gpm_d,gopm_d
  end interface
  interface gcon
     module procedure gcon,gtcon,gcon_d,gtcon_d
  end interface
  interface grect
     module procedure grect0,grect4
  end interface
  interface grange
     module procedure grange1,grange2,gtrange
  end interface
  interface ggappen
     module procedure ggappen,ggoappen
  end interface
interface ggset_status; module procedure ggset_status; end interface
  interface gxax
     module procedure gxax,gtxax
  end interface
  interface gyax
     module procedure gyax,gtyax
  end interface
  interface garr
     module procedure garr
  end interface
contains
real function ggxwtod(t) ! conversions between window and device coordinates
  use pfstatus
  implicit none
  real,intent(IN):: t
  ggxwtod=ax+bx*t;   return
end function ggxwtod
real function ggywtod(t)
  use pfstatus
  implicit none
  real,intent(IN):: t
  ggywtod=ay+by*t; 
end function ggywtod
real function ggxdtow(t); 
  use pfstatus
  implicit none
  real,intent(IN):: t
  ggxdtow=(t-ax)/bx
end function ggxdtow
real function ggydtow(t); 
  use pfstatus
  implicit none
  real,intent(IN):: t
  ggydtow=(t-ay)/by;
end function ggydtow
real function ggyttow(t)
  use pfstatus
  implicit none
  real,intent(IN):: t
  if(yrise==0.)then;ggyttow=t*dchh/by; else;ggyttow=t*dchw/by; endif;
end function ggyttow
real function ggxttow(t);
  use pfstatus
  implicit none
  real,intent(IN):: t
  if(xrun==0.) then;ggxttow=t*dchh/bx; else;ggxttow=t*dchw/bx; endif; 
end function ggxttow
real function ggscaleh(); 
  use pfstatus
  implicit none
  ggscaleh=scaleh;
end function ggscaleh
real function ggscalew(); 
  use pfstatus
  implicit none
  ggscalew=scalew;
end function ggscalew
real function ggdchh();   
  use pfstatus
  implicit none
  ggdchh  =dchh
end function ggdchh

real function wstring(buf,l)
  use pfkinds
  use pfconsts
  use pfstatus
  implicit none
!  real wstring
  character(len=*),intent(IN):: buf 
  integer,intent(IN):: l
  integer,parameter:: mfonts=7
  real,parameter:: fontsc=.002
  logical lfirst
  integer,dimension(32:127,mfonts):: kfwid
  integer i,j,ifont
  save
  data lfirst /.true./
  if (lfirst) then
     include 'pathfont.inc'
     do ifont=1,nfonts
        read(4,40)(kfwid(i,ifont),i=32,127)
     enddo
     close(4)
40   format(16i5)
     lfirst=.false.
  endif
  wstring=0.
  do i=1,l
     j=ichar(buf(i:i))
     if(j < 32 .or. j >= 128)stop 'invalid character in wstring'
     wstring=wstring+kfwid(j,ifonts)
  enddo
  wstring=wstring*fontsc
end function wstring

real function gpink()
  gpink=0.08; return
entry gblue();   gblue  =.667; return
entry gorange(); gorange=.1;   return
entry gyellow(); gyellow=.167; return
entry ggreen();  ggreen =.333; return
entry gpurple(); gpurple=.833; return
entry gred();    gred   =0.  ; return

end function gpink
subroutine hpl2(n,x,y)
  use pfstatus
  implicit none
  integer,intent(IN):: n
  real,dimension(n),intent(IN):: x,y
  integer:: nf
  real,allocatable:: xf(:),yf(:)
  if(n<1)return
  if(kfine < 1 .or. kfine > 5)stop 'In gpl; kfine out of bounds'
  nf=1+(n-1)*kfine; if(n <= 2)nf=n
  allocate(xf(nf),yf(nf))
  if(kfine==1 .or. n<=2)then; xf=x; yf=y
  else;                       call ggfline2(n,nf,x,y,xf,yf)
  endif
  call ggpl2(nf,xf,yf)
  deallocate(xf,yf)
end subroutine hpl2
subroutine htpl2(n,x,y,ggtform)
  use pfstatus
  implicit none
  integer,intent(IN):: n
  real,dimension(n),intent(IN):: x,y
  external ggtform
  integer:: i,nf
  real:: xt,yt
  real,allocatable:: xf(:),yf(:),zf(:)
  if(n<1)return; if(kfine < 1 .or. kfine > 5)stop 'In gpl; kfine out of bounds'
  nf=1+(n-1)*kfine
  allocate(xf(nf),yf(nf),zf(nf))
  if(kfine==1)then;   xf=x; yf=y
  elseif(n <= 2)then; call gglline2(n,nf,x,y,xf,yf)
  else;               call ggfline2(n,nf,x,y,xf,yf)
  endif
  do i=1,nf;xt=xf(i); yt=yf(i); call ggtform(xt,yt,xf(i),yf(i),zf(i)); enddo
  call gopl2(nf,xf,yf,zf)
  deallocate(xf,yf,zf)
end subroutine htpl2

subroutine hpl3(n,x,y,z)
  use pfstatus
  implicit none
  integer,intent(IN):: n
  real,dimension(n),intent(IN):: x,y,z
  integer:: nf
  real,allocatable:: xf(:),yf(:),zf(:)
  if(n<1)return
  if(kfine < 1 .or. kfine > 5)stop 'In gpl; kfine out of bounds'
  nf=1+(n-1)*kfine; if(n <= 2)nf=n
  allocate(xf(nf),yf(nf),zf(nf))
  if(kfine==1 .or. n<=2)then;   xf=x; yf=y; zf=z
  else;                         call ggfline3(n,nf,x,y,z,xf,yf,zf)
  endif
  call gopl2(nf,xf,yf,zf)
  deallocate(xf,yf,zf)
end subroutine hpl3

subroutine hpl2_d(n,x_d,y_d)
  use pfkinds
  implicit none
  integer,              intent(IN):: n
  real(dp),dimension(n),intent(IN):: x_d,y_d
  real,dimension(n):: x,y
  x=x_d; y=y_d
  call hpl2(n,x,y)
end subroutine hpl2_d
subroutine htpl2_d(n,x_d,y_d,ggtform)
  use pfkinds
  implicit none
  integer,              intent(IN):: n
  real(dp),dimension(n),intent(IN):: x_d,y_d
  external ggtform
  real,dimension(n):: x,y
  x=x_d; y=y_d
  call htpl2(n,x,y,ggtform)
end subroutine htpl2_d
subroutine hpl3_d(n,x_d,y_d,z_d)
  use pfkinds
  implicit none
  integer,intent(IN):: n
  real(dp),dimension(n),intent(IN):: x_d,y_d,z_d
  real,dimension(n):: x,y,z
  x=x_d; y=y_d; z=z_d
  call hpl3(n,x,y,z)
end subroutine hpl3_d

subroutine gpm(n,x,y)
  use pfconsts
  use pfstatus
  implicit none
  integer,intent(IN):: n
  real,dimension(n),intent(IN):: x,y
  integer:: i1
  ilpen=impen
  if(n <= 0)return
  do i1=1,n; call gmarker(x(i1),y(i1)); enddo
end subroutine gpm

subroutine gopm(n,x,y,z)
  use pfconsts
  use pfstatus
  implicit none
  integer,intent(in):: n
  real,dimension(n),intent(in):: x,y,z
  integer:: i1
  ilpen=impen
  if(n <= 0)return
  do i1=1,n; if(z(i1) >= 0.)call gmarker(x(i1),y(i1)); enddo
end subroutine gopm

subroutine gpm_d(n,x_d,y_d)
  use pfkinds
  implicit none
  integer,intent(IN):: n
  real(dp),dimension(n),intent(IN):: x_d,y_d
  real,dimension(n):: x,y
  x=x_d; y=y_d
  call gpm(n,x,y)
end subroutine gpm_d

subroutine gopm_d(n,x_d,y_d,z_d)
  use pfkinds
  implicit none
  integer,intent(IN):: n
  real(dp),dimension(n),intent(IN):: x_d,y_d,z_d
  real,dimension(n):: x,y,z
  call gopm(n,x,y,z)
end subroutine gopm_d

subroutine gcon(a)
!  contour using existing contour range and intervals
!  and using present contouring window
  use pfstatus
  implicit none
  real,dimension(:,:),intent(IN):: a
  integer:: j
  real:: x0,x1,y0,y1,xh,yh
  ltform=.false. ! Tells gtcon that NO user-supplied transformation occurs
  call gqnt(j); call gqwn(j,x0,x1,y0,y1)
  call grect(x0,x1,y0,y1)
  call gtcon(a,gguform)
  xh=(x0+x1)/2
  yh=(y0+y1)/2
  call ggscapxy(xh,y0,xh,y1,x0,yh,x1,yh)
  call gsdvxy(5000.,7000.)
  ltform=.true.  ! Revert to the default
end subroutine gcon
subroutine gcon_d(a_d)
!  contour using existing contour range and intervals
!  and using present contouring window
  use pfkinds
  use pfstatus
  implicit none
  real(dp),dimension(:,:),intent(IN):: a_d
  integer:: j
  real:: x0,x1,y0,y1,xh,yh
  ltform=.false. ! Tells gtcon that NO user-supplied transformation occurs
  call gqnt(j); call gqwn(j,x0,x1,y0,y1)
  call grect(x0,x1,y0,y1)
  call gtcon_d(a_d,gguform)
  xh=(x0+x1)/2
  yh=(y0+y1)/2
  call ggscapxy(xh,y0,xh,y1,x0,yh,x1,yh)
  call gsdvxy(5000.,7000.)
  ltform=.true.  ! Revert to the default
end subroutine gcon_d

subroutine gtcon_d(a_d,ggtform)
  use pfkinds
  implicit none
  real(dp),dimension(0:,0:),intent(IN):: a_d
  external ggtform
  integer:: nx,ny
  real,allocatable :: a(:,:)
  nx=size(a_d,1); ny=size(a_d,2)
  allocate (a(0:nx-1,0:ny-1))
  a=a_d
  call gtcon(a,ggtform)
  deallocate(a)
end subroutine gtcon_d

subroutine gtcon(a,ggtform) !  transform and contour
  use pfconsts
  use pfstatus
  implicit none
  real,dimension(0:,0:),intent(IN):: a
  integer:: nx,ny,nxf,nyf,nxm,nym,nxfm,nyfm,iy,ixf
  real:: x0,x1,y0,y1
  integer:: nxy,k,k1,k2
  real:: a1,a2,ad,cmin,cmax

  real,allocatable :: af(:,:),aff(:,:)
  external ggtform
  if(ltform.and..not.lcwn)stop 'In gcon; contour domain window undefined'
  nx=size(a,1); ny=size(a,2)
  print '('' nx,ny='',2i5)',nx,ny
  if(nx<2 .or. ny<2)stop 'In gcon; nx or ny out of bounds'
  if(kfine<1 .or. kfine>5)stop 'In gcon; kfine out of bounds' 
  if(lcwn)then
     call gqcwn(x0,x1,y0,y1)
  else
     call gqwn(itnr,x0,x1,y0,y1)
  endif
  if(.not.ltform.and.lcwn)then
     xw0_old=xw0(itnr);xw1_old=xw1(itnr);yw0_old=yw0(itnr);yw1_old=yw1(itnr)
     call gswn(itnr,x0,x1,y0,y1)
     call gsnt(itnr)
  endif
  if(.not.lcval)then
  call grange(a,a1,a2)
  ncon=0
  ad=a2-a1
  if(ad==0.)then
     print '('' In gcon; no increment defined because range vanishes'')'
     ncon=0; goto 300
  endif
  call gdcade(ad,mcon,dinck,modec)
  ad=a2-a1
  dincj=5*dinck
  k1=nint(a1/dinck+.5)
  k2=nint(a2/dinck-.5)
  cmin=k1*dinck
  cmax=k2*dinck
  do k=k1,k2
     if(mod(k,2)/=0.or.modec==0)then
        ncon=ncon+1
        if(ncon>ncon_max)stop 'In gcon; default number of contours too large'
        cval(ncon)=k*dinck
        if(mod(k,5)==0)then
           if(k>=0)then
              ctyp(ncon)=ipos5
           else
              ctyp(ncon)=ineg5
           endif
        else
           if(k>0)then
              ctyp(ncon)=ipos1
           else
              ctyp(ncon)=ineg1
           endif
        endif
     endif
  enddo
  print '('' in gtcon, actual and intended,ncon,mcon='',2i5)',ncon,mcon
  endif
300 continue

  nxm=nx-1; nym=ny-1
  if(kfine==1 .or. nx == 2 .or. ny == 2)then
     call ggkrcon(a,x0,x1,y0,y1,nx,ny,nx,ny,ggtform)
  else
     if(nx==2)then; nxf=2
     else;          nxf=1+nxm*kfine
     endif
     if(ny==2)then; nyf=2
     else;          nyf=1+nym*kfine
     endif
     nxfm=nxf-1; nyfm=nyf-1
     allocate(af(0:nxfm,0:nym),aff(0:nxfm,0:nyfm))
     if(nx==2)then; af=a
     else; do iy=0,nym; call ggfline1(nx,nxf,a(:,iy),af(:,iy)); enddo
     endif
     if(ny==2)then; aff=af
     else; do ixf=0,nxfm; call ggfline1(ny,nyf,af(ixf,:),aff(ixf,:)); enddo
     endif
     call ggkrcon(aff,x0,x1,y0,y1,nxf,nyf,nxf,nyf,ggtform)
     deallocate(af,aff)
  endif
  if(.not.ltform.and.lcwn)then
     x0=xw0_old;x1=xw1_old;y0=yw0_old;y1=yw1_old
     call gswn(itnr,x0,x1,y0,y1)
     call gsnt(itnr)
  endif
end subroutine gtcon

subroutine grect0
  integer:: j
  real:: x0,x1,y0,y1
  call gqnt(j);call gsnt(j);call gqwn(j,x0,x1,y0,y1);call grect4(x0,x1,y0,y1)
end subroutine grect0

subroutine grect4(x0,x1,y0,y1)
  use pfstatus
  real,intent(IN):: x0,x1,y0,y1
  integer:: j
  real,dimension(0:4):: x,y
  call gginixax
  x(0)=x0;x(1)=x1;x(2)=x1;x(3)=x0;x(4)=x0
  y(0)=y0;y(1)=y0;y(2)=y1;y(3)=y1;y(4)=y0
  j=kfine; kfine=1; call gpl(5,x,y); kfine=j
  call ggfinxax
end subroutine grect4

subroutine grange1(a,a1,a2)!  find min (a1) and max (a2) of array of data a
  implicit none
  real,dimension(:),intent(IN):: a
  real,intent(OUT):: a1,a2
  integer:: nx
  integer:: ix
  real:: ai
  nx=size(a,1)
  a1=a(1);  a2=a(1)
  do ix=1,nx
     ai=a(ix); a1=min(a1,ai); a2=max(a2,ai)
  enddo
end subroutine grange1

subroutine grange2(a,a1,a2)!  find min (a1) and max (a2) of array of data a
  implicit none
  real,dimension(:,:),intent(IN):: a
  real,intent(OUT):: a1,a2
  integer:: nx,ny
  integer:: ix,iy
  real:: ai
  nx=size(a,1); ny=size(a,2)
  a1=a(1,1);  a2=a(1,1)
  do iy=1,ny; do ix=1,nx
     ai=a(ix,iy); a1=min(a1,ai); a2=max(a2,ai)
  enddo; enddo
end subroutine grange2

subroutine gtrange(a,b,c1,c2,d1,d2,ggtform)
  implicit none
  real,dimension(:,:),intent(IN) :: a,b
  real,               intent(OUT):: c1,c2,d1,d2
  integer:: ix,iy, nx,ny
  real::    c,d
  external ggtform
  nx=size(a,1);ny=size(a,2)
  if(nx/=size(b,1).or.ny/=size(b,2))stop 'In grange; dimensions inconsistent'
  call ggtform(a(1,1),b(1,1),c,d)
  c1=c; c2=c; d1=d; d2=d
  do iy=1,ny; do ix=1,nx
     call ggtform(a(ix,iy),b(ix,iy),c,d)
     c1=min(c,c1); c2=max(c,c2); d1=min(d,d1); d2=max(d,d2)
  enddo; enddo
end subroutine gtrange

subroutine ggset_status !=====================================================
use pfkinds
use pfconsts
use pfstatus
call ggindef_a
chh=chh_def
dchh=chh*scaleh; dchw=chh*scalewh
xc2=xc2_def; yc2=yc2_def;
xrun=xrun_def; yrise=yrise_def
xw0=xw0_def; xw1=xw1_def; yw0=yw0_def; yw1=yw1_def
xv0=xv0_def; xv1=xv1_def; yv0=yv0_def; yv1=yv1_def
mtypei=mtypei_def; mpeni=mpeni_def
iframe=0
ifonts=1;
kfine_a=0
ggunit=ggunit_def
imtype=0
impen=0
iltype=0
ilpen=0
itnr=0
iltp=0
imtp=0
itick=0
itpen=0
igpen=0
ilpen_old=0
iltype_old=0
itpen_old=0
jv=0

rpli=1.
rpmi=1.
xcapb=0.; ycapb=0.; xcapl=0.;ycapl=0.; xcapt=0.;ycapt=0.; xcapr=0.;ycapr=0.
zarr=1.

lini=.false. ! Usually redundant (data statement in pfigm1.f90)
lportrait=.false.
lopen=.false.
lgmess=lgmess_def
lstroke=.false.
llabel=.true.
ltform=.true.

fistem=fistem_def
call gdbdiv  !	 default number of boundary line divisions (20)
call gdfine  !	 default curve/contouring refinement factor (=1)
call gdmcon  !	 default max contours	   (=50)
call gdmodec !	 default mode of contour selection (integer multiples)
call gdc15   !	 default contouring line-types for 1's and 5's
call gdcwn   !	 contour window is window of normalized transformation.
call gdcons  !	 default present number of contours (=0)
call gdchh   !	 default character height =1.
call gdpmh   !	 default polymarker height=1.
call gdplh   !	 default line-dash size
call gdpmr   !	 default polymarker representations
call gdplr   !	 default polyline representations
call gddtr   !   default dash-type representations (for polylines)
call gdcval  !   Let contours be chosen from the data according to defaults
call gdpnr
call gdstroke
do k=0,ntnr
   call gdwn(k) !	 default window
   call gdvp(k) !	 default viewport
enddo
call gdsty   !	 default shading type
call gdsdar  !	 default shading dark fraction in broken-shading
call gdspen  !	 default shading pen
call gdsodd  !	 default shading staggering for odd lines.
call gdnt    !   default normalized transformation (0)
pwid=wid_def
lini=.true.
end subroutine ggset_status

subroutine ggscapb(xcb,ycb)
  use pfstatus
                        xcapb=ax+bx*xcb; ycapb=ay+by*ycb; return
entry ggscapt(xct,yct); xcapt=ax+bx*xct; ycapt=ay+by*yct; return
entry ggscapl(xcl,ycl); xcapl=ax+bx*xcl; ycapl=ay+by*ycl; return
entry ggscapr(xcr,ycr); xcapr=ax+bx*xcr; ycapr=ay+by*ycr; return
end subroutine ggscapb

subroutine ggoappen(x,y,z)
  use pfconsts
  use pfstatus
  implicit none
  real,intent(IN):: x,y,z
  jv=jv+1; if(jv > jv_max)stop 'In ggoappen; jv exceeds bound'
  xarr(jv)=x; yarr(jv)=y; zarr(jv)=z;
end  subroutine ggoappen
subroutine ggappen(x,y)
  use pfconsts
  use pfstatus
  implicit none
  real,intent(IN):: x,y
  jv=jv+1; if(jv > jv_max)stop 'In ggappen; jv exceeds bound'
  xarr(jv)=x; yarr(jv)=y; zarr(jv)=1.
end subroutine ggappen

subroutine ggcurve(a,ix,iy,id,ggtform)
!=============================================================================
!  CREATE A COMPLETE NEW CONTOUR, ALTERING KT IN THE PROCESS SO THAT
!  THIS CONTOUR IS NOT DRAWN A SECOND TIME. USE DIRECTION CODE KD TO
!  RECORD THE POSITION OF EACH FRESH GRID SQUARE RELATIVE TO PREVIOUS.
!  KD	    : 1 2 3 4
!  DIRECTION: E N W S
!=============================================================================
  use pfconsts
  use pfstatus
  implicit none
  real,dimension(mnnx,tmy),intent(IN):: a
  integer,intent(in):: ix,iy,id
  integer:: k,jx,jy,jxp,jyp,kd,ktj
  integer,dimension(4) :: kbed7,k1284
  integer,dimension(14):: kdnext
  real::x,y,xc,yc,zc
!  interface
!     subroutine ggtform(x1,y1,x2,y2,z2);
!       	implicit none;
!	real:: x1,y1,x2,y2,z2;
!     end subroutine ggtform
!  end interface
  data  kbed7/11,14,13,7/,k1284/1,2,8,4/
  data kdnext/4,1,1,3,4,0,1,2,0,2,2,3,4,3/
  jx=ix; jy=iy
  kd=id
  do k=1,jv_max
     select case(kd)
     case(1,3)
        x=(jx+px)*sx; y=(jy+py+ggrcon(a(jx,jy),a(jx,jy+1),acon1))*sy
        call ggtform(x,y,xc,yc,zc)
        call ggoappen(xc,yc,zc)
        if(kd==3)jx=jx-1 ! If westward, decrement jx
        if(jx==0.or.jx==tnx)exit
     case(2,4)
        x=(jx+px+ggrcon(a(jx,jy),a(jx+1,jy),acon1))*sx; y=(jy+py)*sy
        call ggtform(x,y,xc,yc,zc)
        call ggappen(xc,yc,zc)
        if(kd==4)jy=jy-1 ! If southward, decrement jy
        if(jy==0.or.jy==tny)exit
     end select
     ktj=kt(jx,jy); if(ktj==0)exit
     if(ktj==6 .or.ktj==9)then
        jxp=jx+1; jyp=jy+1
        ktj=k1284(kd); 
        if(a(jx,jy)+a(jxp,jy)+a(jx,jyp)+a(jxp,jyp)>acon1*4.)ktj=kbed7(kd)
     endif
     kty(jx)=kty(jx)-1
     kt(jx,jy)=kt(jx,jy)-ktj; if(kt(jx,jy)<0)kt(jx,jy)=kt(jx,jy)+15
     kd=kdnext(ktj); if(kd==1)jx=jx+1; if(kd==2)jy=jy+1
  enddo
  call ggoenlin 
  contains
    real function ggrcon(a1,a2,ac); 
      implicit none
      real,intent(IN)::a1,a2,ac
      ggrcon=(ac-a1)/(a2-a1); 
    end function ggrcon
end subroutine ggcurve

subroutine ggkrcon(a,x1,x2,y1,y2,mxt,myt,nnx,nny,ggtform)!  CONTOUR A PATCH
  use pfstatus
  implicit none
  integer,intent(IN):: mxt,myt,nnx,nny
  real,DIMENSION(nnx,nny),intent(IN):: a
  real,intent(IN):: x1,x2,y1,y2
!---------------------------------------------------------------------------
  integer:: icon,tmxm,tmym,ix,iy,lx,ly 
  real:: a1,a2,axy
  external ggtform
!============================================================================
  mnnx=nnx
  tmx=mxt
  tmy=myt
  tmxm=tmx-1
  tmym=tmy-1
  sx=(x2-x1)/tmxm
  sy=(y2-y1)/tmym
  a1=a(1,1)
  a2=a(2,2)
  do ix=1,mxt
     do iy=1,myt
        axy=a(ix,iy)
        a1=min(a1,axy)
        a2=max(a2,axy)
     enddo
  enddo
  do icon=1,ncon
     acon1=cval(icon)
     if(a1 < acon1 .and. acon1 < a2)then
        call gspli(ctyp(icon))
        acon4=acon1*4
        do lx=1,tmxm,32
           tnx=min(33,tmx+1-lx)
           tnxm=tnx-1
           px=x1/sx+lx-2
           do ly=1,tmym,32
              tny=min(33,tmy+1-ly)
              tnym=tny-1
              py=y1/sy+ly-2
              call gglrcon(a(lx,ly),ggtform)
           enddo
        enddo
     endif
  enddo
end subroutine ggkrcon

!============================================================================
!                SUBROUTINE GGLRCON
!  SET UP THE ARRAY KT OF TRAVERSAL CODES. KT RANGES FROM 0 TO 15
!  ACCORDING TO THE LEVELS, RELATIVE TO CONTOUR VALUE ACON1, OF THE
!  VALUES OF A AT THE FOUR CORNERS OF THE GRID-SQUARE IN QUESTION
!
!  DENOTING EDGES OF THE GRID SQUARE "W","S","E","N" AND DENOTING THE
!  CLASS OF CONTOUR SEGMENTS THAT ENTER FROM W, EXIT VIA S, BY "WS" etc.,
!  THEN THE FULL SET OF TRAVERSAL CODES ARE DEFINED AS FOLLOWS:
!  KT		  : 0  1  2  3	4  5	   6	  7
!  EDGES CONNECTED: -- WS SE WE NW NS NW&SE/NE&SW NE
!  KT		  : 15 14 13 12 11 10	   9	  8
!  EDGES CONNECTED: -- SW ES EW WN SN WN&ES/EN&WS EN
!============================================================================
subroutine gglrcon(a,ggtform)
  use pfstatus
  real,dimension(mnnx,tmy),intent(IN):: a
  integer,dimension(0:15):: l139b,l89cd,l26ae,l4567, kdt
  external ggtform
!             0,1,2,3, 4,5,6,7, 8,9,a,b, c,d,e,f
!-------------------------------------------------
  data  l139b/0,1,0,1, 0,0,0,0, 0,1,0,1, 0,0,0,0/&  ! <- W-  (1)
       ,l89cd/0,0,0,0, 0,0,0,0, 1,1,0,0, 1,1,0,0/&  ! <- E-  (3)
       ,l26ae/0,0,1,0, 0,0,1,0, 0,0,1,0, 0,0,1,0/&  ! <- S-  (2)
       ,l4567/0,0,0,0, 1,1,1,1, 0,0,0,0, 0,0,0,0/   ! <- N-  (4)
  data kdt/0,1,1,1,1,1,2,1 , 1,2,1,1,1,1,1,0/
  do ix=1,tnx
     do iy=1,tny; kt(ix,iy)=0; if(a(ix,iy)>acon1)kt(ix,iy)=1; enddo
  enddo
  ktxy=0
  do ix=1,tnxm; ixp=ix+1
     ktys=0
     do iy=1,tnym; iyp=iy+1
        kt(ix,iy)=kt(ix,iy)+2*kt(ixp,iy)+4*kt(ix,iyp)+8*kt(ixp,iyp)
        ktys=ktys+kdt(kt(ix,iy))
     enddo
     kty(ix)=ktys   ! Total number of line segments in column ix of patch
     ktxy=ktxy+ktys ! Total number of line segments in this patch
  enddo
  if(ktxy>0)then
!  contours entering in directions 1 and 3:
      do iy=1,tnym
       if(l139b(kt(   1,iy))==1)call ggcurve(a,  1,iy,1,ggtform)
       if(l89cd(kt(tnxm,iy))==1)call ggcurve(a,tnx,iy,3,ggtform)
      enddo
!  contours entering in directions 2 and 4:
      do ix=1,tnxm
       if(l26ae(kt(ix, 	 1))==1)call ggcurve(a,ix,  1,2,ggtform)
       if(l4567(kt(ix,tnym))==1)call ggcurve(a,ix,tny,4,ggtform)
      enddo
!  interior contours:
      do ix=1,tnx-2
         if(kty(ix)>0)then
            do iy=1,tny-2
               if((kt(ix,iy)==7))call ggcurve(a,ix,iy+1,4,ggtform)
               if((kt(ix,iy)==8))call ggcurve(a,ix+1,iy,3,ggtform)
            enddo
         endif
      enddo
   endif
end subroutine gglrcon

subroutine ggindef_a !========================================================
  use pfstatus
  integer:: i
  real::u
  data u/-999999./; data i/-999999/
  xc2_a=u; yc2_a=u;wid_a=u; hue_a=u; sat_a=u; bri_a=u; chh_a=u
  ifonts_a=i 
  call gdstroke
end subroutine ggindef_a

subroutine gclr ! clear current file-frame and write it to disk
  use pfconsts
  use pfstatus
  if(lopen)then
     if(lstroke)call ggfmt108(ggunit)
     call ggfmt199(ggunit)
     close(ggunit)
     lopen=.false.
     xrun=1.; yrise=0.
     itnr=0
     xc2=xc2_def; yc2=yc2_def
     call ggindef_a
  endif
end subroutine gclr

subroutine ggopen
  use pfconsts
  use pfstatus
  iframe=iframe+1
  if(iframe > iframe_max)stop 'Frame index exceeds maximum'
  if(iframe < 10)then
     write(unit=finame,fmt='(a4,''0'',i1,''.ps'')')fistem,iframe
  elseif(iframe < 100)then
     write(unit=finame,fmt='(a4,i2,''.ps'')')fistem,iframe
  else
     stop 'iframe does not fall within allowed range'
  endif
  open(unit=ggunit,file=finame,access='sequential',status='unknown')
  if(lgmess)print '('' Open graphics output file '',A9)',finame
  lopen=.TRUE.
  if(lportrait)then
     call ggfmt98(ggunit)
  else
     call ggfmt99(ggunit)
  endif
  call ggindef_a
end subroutine ggopen

!=============================================================================
!          SUBROUTINE GGSQUASH
! Remove leading and trailing blanks from a %-terminated string of text
! and return the resulting length, including the %
!=============================================================================
subroutine ggsquash(buf,ltext)
  implicit none
  character*80,intent(inout):: buf
  integer,       intent(out):: ltext
  integer:: i,ib1,ib2,ib3
  ib1=0
  do i=1,80 		         ! loop until
     if(buf(i:i) /= ' ')goto 301 ! non-blank
     ib1=i
  enddo
  ltext=80
  print '('' warning: no percent terminator in ggsquash'')'
  buf(1:1)='%'
  return
301 ib3=ib1
  do i=ib1+1,80		       ! loop until
     if(buf(i:i)=='%')goto 302 ! percent
     ib3=i
  enddo
302 ib2=ib3
  do i=ib3,ib1+1,-1
     ib2=i
     if(buf(i:i) /= ' ')goto 303
  enddo
303 ltext=ib2-ib1
  do i=1,ltext
     buf(i:i)=buf(i+ib1:i+ib1)
  enddo
  buf(ltext+1:ltext+1)='%'
end subroutine ggsquash

subroutine gini
  use pfstatus
             lcolor=.false.; call ggset_status; return
entry ginic; lcolor=.true.;  call ggset_status
end subroutine gini

subroutine gportrait
  use pfstatus
  lportrait=.true.
  return
entry glandscape
  lportrait=.false.
end subroutine gportrait

!-----------------------------------------------------------------------------
!			SUBROUTINE GGDECADE				      
!  Given end points of a line segment, determine appropriate "round number"   
!  increments by which the segment can be nicely divided up into a graticule  
!									      
! --> X1,X2	end points of the line segment				      
! --> AM	an amplifying constant, expressed as its log_10 	      
! <-- IT1	decadal magnitude-range associated with the segment's length  
! <-- RT1	decadal residual					     
! <-- DINC5	coarse increment (5 of the fine increments)		      
! <-- DINC1	fine increment						      
! <-- IX1O5	starting index of coarse increments			      
! <-- IX2O5	ending index of coarse increments			      
! <-- IX1O1	starting index of fine increment			      
! <-- IX2O1	endinf index of fine increment				      
!-----------------------------------------------------------------------------
subroutine ggdecade(x1,x2,am,it1,rt1,dinc5,dinc1,ix1o5,ix2o5,ix1o1,ix2o1)
  implicit none
  real,   intent(IN) :: x1,x2,am
  integer,intent(OUT):: it1,ix1o5,ix2o5,ix1o1,ix2o1
  real,   intent(OUT):: rt1,dinc5,dinc1
  real:: dp,xt1,xt2,at1
  dp=(x2-x1)*1.00001   ! range+bonus
  xt1=x2-dp 	   ! ensure inclusion of end-points
  xt2=x1+dp 	   ! allowing for possible round-off
  at1=alog10(dp)+am
  it1=nint(at1)	   ! find which decadal range applies..
  rt1=at1-it1
  if(rt1 <= 0.)then
     dinc5=10.**it1	   ! ..and whether units..
  else
     dinc5=5*10.**it1	   ! ..or fives are appropriate.
  endif
  dinc1=dinc5/5.	   ! divide main (j) interval into fifths (k)
  ix1o5=nint(xt1/dinc5+.5)
  ix2o5=nint(xt2/dinc5-.5)
  ix1o1=nint(xt1/dinc1+.5)
  ix2o1=nint(xt2/dinc1-.5)
end subroutine ggdecade

subroutine gddtr
  use pfconsts
  use pfstatus
  implicit none
  integer:: i1,i2,i3,i
  ndtlen=1; ndtlen(4)=2; ndtlen(6)=3
  dtlen1=0.; dtlen2=0.
  dtlen1(1,1) =13.
  dtlen1(1,2) =200.
  dtlen1(1,3) =40.
  dtlen1(1,4) =40.; dtlen1(2,4)=200.
  dtlen1(1,5) =120.
  dtlen1(1,6) =40.; dtlen1(2,6)=120.; dtlen1(3,6)=40.
  dtlen1(1,7) =40.
  dtlen1(1,8) =40.
  dtlen1(1,9) =800.
  dtlen1(1,10)=60.

  dtlen2(1,1) =40.
  dtlen2(1,2) =320.
  dtlen2(1,3) =80.
  dtlen2(1,4) =80.; dtlen2(2,4)=320.
  dtlen2(1,5) =200.
  dtlen2(1,6) =80.; dtlen2(2,6)=160.; dtlen2(3,6)=80.
  dtlen2(1,7) =320.
  dtlen2(1,8) =800.
  dtlen2(1,9) =880.
  dtlen2(1,10)=120.
  do i1=11,nltype,10
     i2=min(i1+9,nltype)
     do i3=i1,i2
        i=i3+1-i2
        ndtlen(i3)  =ndtlen(i)
        dtlen1(:,i3)=dtlen1(:,i)
        dtlen2(:,i3)=dtlen2(:,i)
     enddo
  enddo
end subroutine gddtr

subroutine gopl2(n,x,y,z)
  implicit none
  integer,intent(IN):: n
  real,dimension(n),intent(IN):: x,y,z
  logical l1,l2
  integer:: i,n1
  l1=.false.
  do i=1,n
     l2=(z(i) >= 0.)
     if(l1.and..not.l2)call ggpl2(i-n1,x(n1),y(n1))
     if(.not.l1.and.l2)n1=i
     l1=l2
  enddo
  if(l1)call ggpl2(n+1-n1,x(n1),y(n1))
end subroutine gopl2

subroutine ggpl2(n,xc,yc)
use pfconsts
use pfstatus
implicit none
integer,          intent(IN):: n
real,dimension(n),intent(IN):: xc,yc
integer:: i100,i100m,j100,i,j,ic1,ic2,id,iseg,nc,k
real:: sc1,sc2,xc1,yc1,wd,w2,snow,sdend,wid,hue,sat,bri
real,dimension(101):: xd,yd
if(n<=1)return
if(.not.lopen)call ggopen
if(iltype<0 .or. iltype>nltype)stop 'In ggpl2; iltype out of bounds'

wid=pwid(ilpen)
if(wid_a /= wid)then
   wid_a=wid	 
   call ggfmt107(ggunit,wid)
endif
if(lcolor)then
   hue=phue(ilpen); sat=psat(ilpen); bri=pbri(ilpen)
   if(hue_a /= hue .or. sat_a /= sat .or. bri /= bri)then
      hue_a = hue; sat_a = sat; bri_a = bri
      call ggfmt127(ggunit,hue,sat,bri)
   endif
endif
if(iltype==0)then
   do i100=1,n-1,100
      i100m=i100-1
      j100=min(i100+100,n)
      nc=j100-i100m

      do i=1,nc
	 j=i100m+i
	 xd(i)=ax+bx*xc(j)
	 yd(i)=ay+by*yc(j)
      enddo
      call gghpl(ggunit,nc,xd,yd)
   enddo
   xc2=xd(nc); xc2_a=xc2
   yc2=yd(nc); yc2_a=yc2
elseif(iltype > 0 .and. iltype <= nltype)then
   do i100=1,n-1,100
      j100=min(i100+100,n)
      ic1=i100
      ic2=ic1+1
      sc1=0.
      xc1=ax+bx*xc(ic1)
      xc2=ax+bx*xc(ic2)
      yc1=ay+by*yc(ic1)
      yc2=ay+by*yc(ic2)
      wd=sqrt((xc2-xc1)**2+(yc2-yc1)**2)
      sc2=sc1+wd
      wd=1.e-6+wd
      snow=0.
loop_j: do j=1,9999
	 do iseg=1,ndtlen(iltype)
! begin a new dash-segment:
            id=1
            w2=(snow-sc1)/wd
            xd(id)=xc1+w2*(xc2-xc1)
            yd(id)=yc1+w2*(yc2-yc1)
            sdend=snow+rpli*dtlen1(iseg,iltype)
401         if(sdend < sc2)then
               id=id+1
               w2=(sdend-sc1)/wd
               xd(id)=xc1+w2*(xc2-xc1)
               yd(id)=yc1+w2*(yc2-yc1)
               call gghpl(ggunit,id,xd,yd)
               snow=snow+rpli*dtlen2(iseg,iltype)
               do k=1,9999
                  if(snow <= sc2)exit	
                  if(ic2 == j100)exit loop_j
                  ic1=ic2
                  xc1=xc2
                  yc1=yc2
                  sc1=sc2
                  ic2=ic1+1
                  xc2=ax+bx*xc(ic2)
                  yc2=ay+by*yc(ic2)
                  wd=sqrt((xc2-xc1)**2+(yc2-yc1)**2)
                  sc2=sc1+wd
                  wd=1.e-6+wd
               enddo
            else
               id=id+1
               xd(id)=xc2
               yd(id)=yc2
               if(ic2 == j100)then
                  call gghpl(ggunit,id,xd,yd)
                  exit loop_j
               endif
               ic1=ic2
               xc1=xc2
               yc1=yc2
               sc1=sc2
               ic2=ic1+1
               xc2=ax+bx*xc(ic2)
               yc2=ay+by*yc(ic2)
               wd=sqrt((xc2-xc1)**2+(yc2-yc1)**2)
               sc2=sc1+wd
               wd=1.e-6+wd
               goto 401
	    endif
	 enddo
      enddo loop_j
   enddo
endif
lstroke=.true. ! Signifies that "stroke" is now needed in .ps file
end subroutine ggpl2
subroutine gghpl(ggunit,n,x,y)
implicit none
integer,intent(IN):: ggunit,n
real,dimension(n),intent(IN):: x,y
real,dimension(4):: xp,yp
integer:: i,i1,i2,i3,i4
xp(1)=x(1)
yp(1)=y(1)
if(n==2)then  ! this case is so frequent it is worth special code
   call ggfmt106(ggunit,x(1),y(1),x(2),y(2))
   call ggfmt108(ggunit)
else
   call ggfmt100(ggunit,xp(1),yp(1))
   do i1=2,n,4
      i2=min(n,i1+3)
      i3=i2+1-i1
      do i=1,i3
         i4=i1+i-1
         xp(i)=x(i4)
         yp(i)=y(i4)
      enddo
      select case (i3)
         case(1)
            call ggfmt101(ggunit,xp(1),yp(1))  
         case(2)
            call ggfmt102(ggunit,xp(1),yp(1),xp(2),yp(2))
         case(3)
            call ggfmt103(ggunit,xp(1),yp(1),xp(2),yp(2),xp(3),yp(3))
         case(4)
            call ggfmt104(ggunit,xp(1),yp(1),xp(2),yp(2),&
                 xp(3),yp(3),xp(4),yp(4))
      end select
   enddo
   call ggfmt108(ggunit)
endif
end subroutine gghpl

!!!!!!!!!!! Set and query routine with one integer argument:
subroutine gsfine(j)
use pfconsts; use pfstatus
integer:: j
   if(j<1 .or. j>5)stop 'in gsfine; j out of allowed range, [1,5]'
   kfine=j;return;               entry gqfine(j); j=kfine;             return
entry gsmodec(j)
   if(j<0 .or. j>2)stop 'in gsmodec;contouring mode index out of range [0,2]'
   modec=j; return;              entry gqmodec(j); j=modec;            return
entry gsmcon(j); mcon=j; return; entry gqmcon(j);  j=mcon;             return
entry gslabel(j);llabel=j/=0;return;entry gqlabel(j);j=0;if(llabel)j=1;return
entry gstick(j);itick=j; return; entry gqtick(j); j=itick;             return
entry gsnt(j); if(j < 0 .or. j > ntnr)stop 'In gsnt; j out of range'
  itnr=j
  bx=(xv1(itnr)-xv0(itnr))/(xw1(itnr)-xw0(itnr))
  by=(yv1(itnr)-yv0(itnr))/(yw1(itnr)-yw0(itnr))
  ax=xv0(itnr)-bx*xw0(itnr); ay=yv0(itnr)-by*yw0(itnr)
  bx=dy*bx; by=dy*by; ax=dy*ax; ay=dy*ay;                              return
entry gqnt(j); j=itnr;                                                 return
entry gsmess(j); lgmess=.false.; if(j == 0)lgmess=.true.;              return  
entry gqmess(j); j=0;            if(lgmess)j=1;                        return
entry gstpen(j); itpen=j; return; entry gqtpen(j); j=itpen;            return
entry gsgpen(j); igpen=j; return; entry gqgpen(j); j=igpen;            return
end subroutine gsfine

!!!!!!!!!!!!!!!! Set and Query routines with several integer arguments:
subroutine gsc15(jpos1,jneg1,jpos5,jneg5) ! contour types at 1's and 5's
use pfconsts; use pfstatus
integer:: jpos1,jneg1,jpos5,jneg5
   if(jpos1<0 .or. jpos1>nltp)stop 'In gsc15; jpos1 out of bounds'
   if(jneg1<0 .or. jneg1>nltp)stop 'In gsc15; jneg1 out of bounds'
   if(jpos5<0 .or. jpos5>nltp)stop 'In gsc15; jpos5 out of bounds'
   if(jneg5<0 .or. jneg5>nltp)stop 'In gsc15; jneg5 out of bounds'
   ipos1=jpos1; ineg1=jneg1; ipos5=jpos5; ineg5=jneg5; return
entry gqc15(jpos1,jneg1,jpos5,jneg5)  ! query contour types
   jpos1=ipos1; jneg1=ineg1; jpos5=ipos5; jneg5=ineg5; return
end subroutine gsc15
subroutine gsdtr(j,nglen,glen1,glen2)
!  reset line-dash-type j to a user-defined line-type
  use pfconsts
  use pfstatus
  implicit none
  integer,intent(IN)             :: j
  integer,intent(INOUT)          :: nglen
  real,dimension(3),intent(INOUT):: glen1,glen2
  integer:: i
  if(nglen<1.or.nglen>3)stop 'In gsdtr; nglen out of range'
  if(j<1.or.j>nltype)     stop 'In gsdtr; j out of range'
  ndtlen(j)=nglen
  if(lgmess)then
     print '('' linetype'',i2,'' redefined'')',j
     print '('' new segment parameters are:'')'
     do i=1,nglen; print '(2(1x,f8.0))',glen1(i),glen2(i); enddo
  endif
  do i=1,3
     dtlen1(i,j)=glen1(i)
     dtlen2(i,j)=glen2(i)
  enddo
  return
entry gqdtr(j,nglen,glen1,glen2)
  if(j<1.or.j>nltype)     stop 'In gqdtr; j out of range'
  nglen=ndtlen(j)
  do i=1,nglen
     glen1(i)=dtlen1(i,j)
     glen2(i)=dtlen2(i,j)
  enddo
end subroutine gsdtr
subroutine gspmr(j,jtype,jpen)
  use pfstatus
  implicit none
  integer,intent(IN):: j
  integer:: jtype,jpen
  if(j==0.and.lini)stop 'In gspmr; j=0 forbidden except at initialization'
  if(jtype<0 .or. jtype>nmtype)stop 'In gspmr; jtype out of range'
  if(jpen<0  .or. jpen >npen  )stop 'In gspmr; jpen  out of range'
  mtypei(j)=jtype; mpeni(j)=jpen; return
entry gqpmr(j,jtype,jpen); jtype=mtypei(j); jpen=mpeni(j); return
entry gsplr(j,jtype,jpen); 
  if(j==0.and.lini)stop 'In gsplr; j=0 forbidden except at initialization'
  if(jtype<0 .or. jtype>nltype)stop 'In gsplr; jtype out of range'
  if(jpen<0  .or. jpen >npen  )stop 'In gsplr; jpen  out of range'
  ltypei(j)=jtype; lpeni(j)=jpen; return 
entry gqplr(j,jtype,jpen); jtype=ltypei(j); jpen=lpeni(j); return
end subroutine gspmr

!!!!!!!!!!! Set and query various windows:
subroutine gswn(j,x0,x1,y0,y1)
  use pfstatus
  implicit none
  integer,intent(IN):: j
  real:: x0,x1,y0,y1
                           xw0(j)=x0; xw1(j)=x1; yw0(j)=y0; yw1(j)=y1; return
entry gqwn(j,x0,x1,y0,y1); x0=xw0(j); x1=xw1(j); y0=yw0(j); y1=yw1(j); return
entry gsvp(j,x0,x1,y0,y1); xv0(j)=x0; xv1(j)=x1; yv0(j)=y0; yv1(j)=y1; return
entry gqvp(j,x0,x1,y0,y1); x0=xv0(j); x1=xv1(j); y0=yv0(j); y1=yv1(j); return
end subroutine gswn

subroutine gspmh(a)  ! set polymarker height (=size)
  use pfconsts
  use pfstatus
  implicit none
  real:: a
                rpmi=a; return; entry gqpmh(a); a=rpmi; return
entry gsplh(a); rpli=a; return; entry gqplh(a); a=rpli; return
entry gschh(a); chh=a; dchh=chh*scaleh; dchw=chh*scalewh;return
entry gqchh(a); a=chh; return 
end subroutine gspmh
!!!!!!!!!!!!!!!! <== End of Set and query group of routines

subroutine gline(x1,y1,x2,y2)
  use pfstatus
  implicit none
  real,intent(IN):: x1,y1,x2,y2
  if(jv /= 0)call ggenlin
  call ggappen(x1,y1); call ggappen(x2,y2); call ggenlin
end subroutine gline

subroutine ggdump
  use pfconsts
  use pfstatus
  print '('' integer constants:'')'
  print '('' npen          = '',i8)',npen
  print '('' ntnr          = '',i8)',ntnr
  print '('' nltp          = '',i8)',nltp
  print '('' nmtp          = '',i8)',nmtp
  print '('' mtypei_def    = '',i8)',mtypei_def
  print '('' mpeni_def     = '',i8)',mpeni_def
  print '('' ltypei_def    = '',i8)',ltypei_def
  print '('' lpeni_def     = '',i8)',lpeni_def
  print '('' ggunit_def    = '',i8)',ggunit_def
  print '('' nfonts        = '',i8)',nfonts
  print '('' jv_max        = '',i8)',jv_max
  read(*,*)
  print '('' real constants:'')'
  print '('' dy            = '',f18.4)',dy
  print '('' c30           = '',f18.4)',c30
  print '('' s30           = '',f18.4)',s30
  print '('' scaleh        = '',f18.4)',scaleh
  print '('' scalew        = '',f18.4)',scalew
  print '('' scalewh       = '',f18.4)',scalewh
  print '('' chh_def       = '',f18.4)',chh_def
  print '('' xc2_def       = '',f18.4)',xc2_def
  print '('' yc2_def       = '',f18.4)',yc2_def
  print '('' xw0_def       = '',f18.4)',xw0_def
  print '('' xw1_def       = '',f18.4)',xw1_def
  print '('' yw0_def       = '',f18.4)',yw0_def
  print '('' yw1_def       = '',f18.4)',yw1_def
  print '('' xv0_def       = '',f18.4)',xv0_def
  print '('' xv1_def       = '',f18.4)',xv1_def
  print '('' yv0_def       = '',f18.4)',yv0_def
  print '('' yv1_def       = '',f18.4)',yv1_def
  print '('' wid_def       = '',f18.4)',wid_def
  print '('' cfudge        = '',f18.4)',cfudge

  print '('' fistem_def        = '',a4)',fistem_def
  read(*,*)
  print '('' Status variables'')'

  print '('' ax            = '',f18.4)',ax
  print '('' ay            = '',f18.4)',ay
  print '('' bx            = '',f18.4)',bx
  print '('' by            = '',f18.4)',by
  print '('' xc2           = '',f18.4)',xc2
  print '('' yc2           = '',f18.4)',yc2
  print '('' xc2_a         = '',f18.4)',xc2_a
  print '('' yc2_a         = '',f18.4)',yc2_a
  print '('' chh           = '',f18.4)',chh
  print '('' dchh          = '',f18.4)',dchh
  print '('' dchw          = '',f18.4)',dchw
  print '('' chh_a         = '',f18.4)',chh_a
  print '('' xrun          = '',f18.4)',xrun
  print '('' yrise         = '',f18.4)',yrise
  print '('' rpli          = '',f18.4)',rpli
  print '('' rpmi          = '',f18.4)',rpmi
  read(*,*)
  print '('' windows:'')'
  do i=0,ntnr
     write(6,64)i,xw0(i),xw1(i),yw0(i),yw1(i)
  enddo
  read(*,*)
  print '('' viewports:'')'
  do i=0,ntnr
     write(6,64)i,xv0(i),xv1(i),yv0(i),yv1(i)
  enddo
  read(*,*)
  print '('' marker types and pens'')'
  do i=0,nmtp
     write(6,60)i,mtypei(i),mpeni(i)
  enddo
  read(*,*)
  print '('' line types and pens'')'
  do i=0,nltp
     write(6,60)i,ltypei(i),lpeni(i)
  enddo
  read(*,*)
  print '('' pen widths,hues,saturations,brightnesses'')'
  do i=0,npen
     write(6,64)i, pwid(i),phue(i),psat(i),pbri(i)
  enddo
  read(*,*)
  print '('' integer scalar status variables:'')'
  print '('' jv            = '',i8)',jv
  print '('' iframe        = '',i8)',iframe
  print '('' ggunit        = '',i8)',ggunit
  print '('' ifonts        = '',i8)',ifonts
  print '('' ifonts_a      = '',i8)',ifonts_a
  print '('' imtype        = '',i8)',imtype
  print '('' impen         = '',i8)',impen
  print '('' iltype        = '',i8)',iltype
  print '('' ilpen         = '',i8)',ilpen
  print '('' itnr          = '',i8)',itnr
  print '('' iltp          = '',i8)',iltp
  print '('' imtp          = '',i8)',imtp
  print '('' nbdiv         = '',i8)',nbdiv
  print '('' nbdivp        = '',i8)',nbdivp
  print '('' itick         = '',i8)',itick
  print '('' ilpen         = '',i8)',ilpen
  print '('' itpen         = '',i8)',itpen
  print '('' igpen         = '',i8)',igpen
  read(*,*)
  print '('' character variables:'')'
  print '('' fistem        = '',a4)',fistem
  print '('' finame        = '',a9)',finame
  read(*,*)
  print '('' logical variables:'')'
  print '('' lini          = '',l1)',lini
  
  read(*,*)
64 format(1x,i3,4(1x,e12.6))
60 format(1x,i3,4x,4(1x,i8))
end subroutine ggdump

subroutine gspnr(ipen,wid,hue,sat,bri)! Pen representation
use pfconsts
use pfstatus
implicit none
integer,intent(IN):: ipen
real,   intent(IN):: wid,hue,sat,bri
if(ipen < 0 .or. ipen > npen)stop 'In gspnr; ipen out of range'
if(lini .and. ipen == 0)stop 'In gspnr; ipen=0 outside initialization step'
if(wid < 0.)stop 'In gspnr; wid out of range'
if(hue < 0. .or. hue > 1.)stop 'In gspnr; hue out of range'
if(sat < 0. .or. sat > 1.)stop 'In gspnr; sat out of range'
if(bri < 0. .or. bri > 1.)stop 'In gspnr; bri out of range'
pwid(ipen)=wid; phue(ipen)=hue; psat(ipen)=sat; pbri(ipen)=bri
end subroutine gspnr
 
subroutine gdpnr
use pfconsts
use pfstatus
implicit none
integer:: ipen
if(.not.lini)call gspnr(0, wid_def,0.,0.,0.)
if(lcolor)then
   call gspnr(1,  wid_def,gred(),   1.,1.)
   call gspnr(2,  wid_def,gblue(),  1.,1.)
   call gspnr(3,  wid_def,ggreen(), 1.,1.)
   call gspnr(4,  wid_def,gpurple(),1.,1.)
   call gspnr(5,  wid_def,gorange(),1.,1.)
   call gspnr(6,  wid_def,gblue(),  .5,1.)
   call gspnr(7,  wid_def,gred(),   .5,1.)
   call gspnr(8,  wid_def,ggreen(), .5,1.)
   call gspnr(9,  wid_def,gyellow(),1.,1.)
   call gspnr(10, wid_def,gpurple(),.5,1.)
   call gspnr(11, wid_def,gyellow(),1.,.5)
   call gspnr(12, wid_def,gred(),   1.,.5)
   call gspnr(13, wid_def,gblue(),  1.,.5)
   call gspnr(14, wid_def,ggreen(), 1.,.5)
   call gspnr(15, wid_def,gpurple(),1.,.5)
   call gspnr(16, wid_def,gorange(),1.,.5)
   call gspnr(17, wid_def,gyellow(),.5,1.)
   call gspnr(18, wid_def,0.,0.,.5)
   do ipen=19,npen; call gspnr(ipen, wid_def,0.,0.,0.); enddo
else
   do ipen=1,npen; call gspnr(ipen,wid_def,0.,0.,0.); enddo
endif
end subroutine gdpnr
subroutine gdplr
use pfconsts
use pfstatus
implicit none
integer i,j
if(.not.lini)call gsplr(0,0,0)
if(lcolor)then
   j=min(npen,nltp)
   do i=1,j
      call gsplr(i,0,i)
   enddo
   do i=j+1,nltp
      call gsplr(i,0,0)
   enddo
else
   do i=1,nltp; call gsplr(i,0,0); enddo
   call gsplr( 1,1,0) !------------------------------------- 
   call gsplr( 2,2,0) !-----   -----   -----   -----   ----- 
   call gsplr( 3,3,0) !- - - - - - - - - - - - - - - - - - - 
   call gsplr( 4,4,0) !- ---   - ---   - ---   - ---   - --- 
   call gsplr( 5,5,0) !---  ---	---  ---  ---  ---  --- 
   call gsplr( 6,6,0) !- --- - - --- - - --- - - --- - - ---
   call gsplr(11,1,5)  ! used for heavier contours > 0.	
   call gsplr(13,3,5)  ! used for heavier contours < 0.
endif
end subroutine gdplr

subroutine gdpmr
  use pfconsts
  integer:: i
  call gspmr(1,1,0)   !	.  .  . 	      
  call gspmr(2,2,0)   !	*  *  * 
  call gspmr(3,3,0)
  call gspmr(4,4,0)   !	o  o  o 		 
  call gspmr(5,5,0)   !	x  x  x 		 
  call gspmr(6,6,0)   !	  triangles		 
  call gspmr(7,7,0)   !	  diamonds		 
  call gspmr(8,8,0)   !	  squares		 
  call gspmr(9,9,0)   !	  stars of david
  do i=10,nmtp
     call gspmr(i,2,0)  !	+  +  + 		
  enddo
end subroutine gdpmr	

subroutine ggdefaults
integer jtnr
entry gdmess;     call gsmess (1); return ! Message switch (On)
entry gdlabel;    call gslabel(1); return ! Graph labels (On)
entry gdtick;     call gstick (0); return ! Graph axis tick style (Centered)
entry gdnt;       call gsnt   (0); return ! Normalized transformation
entry gdchh;      call gschh (1.); return ! Character height
entry gdplh;      call gsplh (1.); return ! Polyline dash height 
entry gdpmh;      call gspmh (1.); return ! Polymarker height
entry gdwn(jtnr); call gswn(jtnr,0.,1.,0.,1.); return! Window
entry gdvp(jtnr); call gsvp(jtnr,0.,1.,0.,1.); return! Viewport 
entry gdfine;     call gsfine (1); return ! Curve refinement factor
entry gdwnxy;     call gswnxy(0.,0.); return
entry gdtpen;     call gstpen (0); return ! Text pen
entry gdgpen;     call gsgpen (0); return ! Graph frame pen
entry gdmodec;    call gsmodec(0); return ! Contour mode
entry gdc15;      call gsc15(1,3,11,13); return ! Contour types at 1s, 5s.
entry gdmcon;     call gsmcon(50); return
entry gdtang;     call gstang(1.,0.); return ! Text angle
entry gdbdiv;     call gsbdiv(20); return ! Boundary divisions
entry gdfont;     call gsfont (1); return ! Font (Helvetica) 
end subroutine ggdefaults	

subroutine gtcaption(buf)
  use pfstatus
  character(len=*) buf
  call gsdvxy(xcapt,ycapt); call greltx(0.,1.);
  call gqdvxy(xcapt,ycapt); call gcaption(buf); return
entry gbcaption(buf)
  call gsdvxy(xcapb,ycapb); call greltx(0.,-1.5)
  call gqdvxy(xcapb,ycapb); call gcaption(buf); return
entry glcaption(buf)
  call gsdvxy(xcapl,ycapl); call gstang(0.,1.)
  call greltx(0.,.5);       call gqdvxy(xcapl,ycapl)
  call gcaption(buf);       call gdtang;        return
entry grcaption(buf)
  call gsdvxy(xcapr,ycapr); call gstang(0.,1.)
  call greltx(0.,-1.5);     call gqdvxy(xcapr,ycapr)
  call gcaption(buf);       call gdtang
end subroutine gtcaption

subroutine gscwn(x0,x1,y0,y1)
  use pfstatus
  real:: x0,x1,y0,y1
             lcwn=.true.; x0cwn=x0; x1cwn=x1; y0cwn=y0; y1cwn=y1; return
entry gqcwn(x0,x1,y0,y1); x0=x0cwn; x1=x1cwn; y0=y0cwn; y1=y1cwn; return 
entry gdcwn; lcwn=.false.
end subroutine gscwn

subroutine gsfile(cccc)
use pfstatus
implicit none
character*4,intent(IN):: cccc
call gclr
iframe=0
fistem=cccc
if(lgmess)print '('' New output files have name stem: '',a4)',fistem
end subroutine gsfile

subroutine ggenlin
  use pfstatus
  call ggpl2(jv,xarr,yarr)
  jv=0
end subroutine ggenlin
subroutine ggoenlin
  use pfstatus
  call gopl2(jv,xarr,yarr,zarr)
  zarr(1:jv)=1.
  jv=0
end subroutine ggoenlin

subroutine ggupdate
  use pfconsts
  use pfstatus
  implicit none
  real:: wid,hue,sat,bri
  if(.not.lopen)call ggopen
  if(lstroke)call ggfmt108(ggunit)
  if(chh /= chh_a .or. ifonts /= ifonts_a)then
     chh_a=chh
     ifonts_a=ifonts
     call ggfmt112(chh*cfudge)
  endif
  if(wid_a /= pwid(impen))then
     wid =pwid(impen)
     wid_a =wid
     call ggfmt107(ggunit,wid)
  endif
  if(lcolor)then
     hue=phue(ilpen); sat=psat(ilpen); bri=pbri(ilpen)
     if(hue_a /= hue .or. sat_a /= sat .or. bri_a /= bri)then
        hue_a=hue; sat_a=sat; bri_a=bri
        call ggfmt127(ggunit,hue,sat,bri)
     endif
  endif
return
entry ggupdatex
  xc2_a=xc2; yc2_a=yc2
  call ggfmt100(ggunit,xc2,yc2)
end subroutine ggupdate

subroutine gdstroke; use pfstatus; lstroke=.false.; end subroutine gdstroke

subroutine garrow(x1,y1,x2,y2) ! draw arrow from (x1,y1) to (x2,y2) in window
  use pfconsts
  use pfstatus
  implicit none
  real,intent(IN):: x1,y1,x2,y2
  real:: xc1,yc1,xd,yd,sd,xc3,yc3,xc4,yc4
  ilpen_old=ilpen
  ilpen=itpen
  call ggupdate
  xc1=ax+bx*x1
  yc1=ay+by*y1
  xc2=ax+bx*x2
  yc2=ay+by*y2
  xd=xc2-xc1
  yd=yc2-yc1
  sd=sqrt(xd*xd+yd*yd)
  xd=xd*chh*scaleh/sd
  yd=yd*chh*scaleh/sd
  xc3=xc2-c30*xd-s30*yd
  yc3=yc2+s30*xd-c30*yd
  xc4=xc2-c30*xd+s30*yd
  yc4=yc2-s30*xd-c30*yd
  call ggfmt106(ggunit,xc1,yc1,xc2,yc2)
  call ggfmt108(ggunit)
  call ggfmt105(ggunit,xc3,yc3,xc2,yc2,xc4,yc4,xc2,yc2)
  call ggfmt108(ggunit)
  ilpen=ilpen_old
end subroutine garrow

subroutine gtx(xx,yy,buf)	 ! write text starting at (xx,yy)
  use pfconsts
  use pfstatus
  implicit none
  real,intent(IN):: xx,yy
  character*80,intent(IN):: buf
  integer:: l
  real:: wstr
  ilpen=itpen
  call ggupdate
  xc2=ax+bx*xx
  yc2=ay+by*yy
  call ggupdatex
  l=len_trim(buf)		 ! count up the characters of text
  call ggfmt109(ggunit,buf,l)
  wstr=dchw*wstring(buf,l)
  xc2=xc2+wstr*xrun
  yc2=yc2+wstr*yrise
end subroutine gtx

subroutine gsfont(ifont)
  use pfconsts
  use pfstatus
  implicit none
  integer:: ifont
  if(ifont < 1 .or. ifont > nfonts)stop 'In gsfont; ifont out of bounds'
  ifonts=ifont; return
  entry gqfont(ifont); ifont=ifonts;
end subroutine gsfont

subroutine gmarker(x,y)
  use pfconsts
  use pfstatus
  implicit none
  real,intent(IN):: x,y
  real:: xc,yc,xcp,xcm,ycp,ycm,xcpc30,xcmc30,xcps30,xcms30,&
       ycpc30,ycmc30,ycps30,ycms30,rpm,rc30,rs30
  rpm=rpmi*45.
  rc30=rpm*c30; rs30=rpm*s30
  ilpen=impen
  call ggupdate
  xc=ax+bx*x
  yc=ay+by*y
  if(imtype==1)then
     xcp=xc+1.
     xcm=xc-1.
     ycp=yc+1.
     ycm=yc-1.
  else
     xcp=xc+rpm
     xcm=xc-rpm
     ycp=yc+rpm
     ycm=yc-rpm
  endif
  if(imtype==0 .or. imtype==4 .or. imtype==6 .or. imtype>=9)then
     xcpc30=xc+rc30
     xcmc30=xc-rc30
     xcps30=xc+rs30
     xcms30=xc-rs30
     ycpc30=yc+rc30
     ycmc30=yc-rc30
     ycps30=yc+rs30
     ycms30=yc-rs30
  endif
  select case(imtype)
  case(0) ! "+" and "X"
     call ggfmt105(ggunit,xc,ycp, xc,ycm, xcp,yc, xcm,yc)
     call ggfmt105(ggunit,xcp,ycm, xcm,ycp, xcm,ycm, xcp,ycp)
     xc2_a=xcp; yc2_a=ycp
  case(1) ! Tiny square
     call ggfmt100(ggunit,xcp,ycp)
     call ggfmt104(ggunit,xcp,ycm, xcm,ycm, xcm,ycp, xcp,ycp)
     xc2_a=xcp; yc2_a=ycp
  case(2) ! "+"
     call ggfmt105(ggunit,xc,ycp, xc,ycm, xcp,yc, xcm,yc)
     xc2_a=xcm; yc2_a=yc
  case(3) ! "X" and "+"
     call ggfmt105(ggunit,xcp,ycp, xcm,ycm, xcp,ycm, xcm,ycp)
     call ggfmt105(ggunit,xc,ycp, xc,ycm, xcp,yc, xcm,yc)
     xc2_a=xcm; yc2_a=yc
  case(4)!   dodecagon
     call ggfmt100(ggunit,xcp,yc)
     call ggfmt104(ggunit,xcpc30,ycps30,xcps30,ycpc30,xc,ycp,xcms30,ycpc30)
     call ggfmt104(ggunit,xcmc30,ycps30,xcm,yc,xcmc30,ycms30,xcms30,ycmc30)
     call ggfmt104(ggunit,xc,ycm, xcps30,ycmc30, xcpc30,ycms30, xcp,yc)
     xc2_a=xcp; yc2_a=yc
  case(5) ! "X"
     call ggfmt105(ggunit,xcp,ycp, xcm,ycm, xcp,ycm, xcm,ycp)
     xc2_a=xcm; yc2_a=ycp
  case(6) ! Triangle
     call ggfmt100(ggunit,xcmc30,ycms30)
     call ggfmt103(ggunit,xcpc30,ycms30, xc,ycp, xcmc30,ycms30)
     xc2_a=xcmc30; yc2_a=ycms30
  case(7) ! Diamond
     call ggfmt100(ggunit,xcm,yc)
     call ggfmt104(ggunit,xc,ycm, xcp,yc, xc,ycp, xcm,yc)
     xc2_a=xcm; yc2_a=yc
  case(8) ! Square
     call ggfmt100(ggunit,xcm,ycm)
     call ggfmt104(ggunit,xcp,ycm, xcp,ycp, xcm,ycp, xcm,ycm)
     xc2_a=xcm; yc2_a=ycm
  case(9) ! Star of David
     call ggfmt100(ggunit,xcmc30,ycms30)
     call ggfmt103(ggunit,xcpc30,ycms30, xc,ycp, xcmc30,ycms30)
     call ggfmt100(ggunit,xcpc30,ycps30)
     call ggfmt103(ggunit,xcmc30,ycps30, xc,ycm, xcpc30,ycps30)
     xc2_a=xcpc30; yc2_a=ycps30
  case(10) ! Inverted triangle
     call ggfmt100(ggunit,xcpc30,ycps30)
     call ggfmt103(ggunit,xcmc30,ycps30, xc,ycm, xcpc30,ycps30)
     xc2_a=xcpc30; yc2_a=ycps30
  case default
     stop 'In gmarker; imtype out of bounds'
  end select
  lstroke=.true.
  xc2=xc; yc2=yc
end subroutine gmarker
  

subroutine gginixax
  use pfstatus
  implicit none
  if(.not.lopen)call ggopen
  dchh=chh*scaleh; dchw=chh*scalewh
  ilpen_old=ilpen;   ilpen=igpen; 
  iltype_old=iltype; iltype=0
  itpen_old=itpen;   itpen=igpen
  call gdtang	    ! ensure default text orientation
return
entry ggfinxax;  ilpen=ilpen_old; iltype=iltype_old; itpen=itpen_old
end subroutine gginixax

subroutine gglline2(n,nf,x,y,xf,yf)
  use pfstatus
  implicit none
  integer,intent(IN):: n,nf
  real,dimension(0:n-1) :: x,y
  real,dimension(0:nf-1):: xf,yf
  integer:: jk,nm,nmm,nmmkfine
  real:: w,wc
  nm=n-1; nmm=n-2; nmmkfine=nmm*kfine
  xf(0:nf-1:kfine)=x; yf(0:nf-1:kfine)=y
  do jk=1,kfine-1
     w=real(jk)/kfine; wc=1.-w
     xf(jk:jk+nmmkfine:kfine)=wc*x(0:nmm)+w*x(1:nm)
     yf(jk:jk+nmmkfine:kfine)=wc*y(0:nmm)+w*y(1:nm)
  enddo
end subroutine gglline2
subroutine gglline3(n,nf,x,y,z,xf,yf,zf)
  use pfstatus
  implicit none
  integer,intent(IN):: n,nf
  real,dimension(0:n-1) :: x,y,z
  real,dimension(0:nf-1):: xf,yf,zf
  integer:: jk,nm,nmm,nmmkfine
  real:: w,wc
  nm=n-1; nmm=n-2; nmmkfine=nmm*kfine
  xf(0:nf-1:kfine)=x; yf(0:nf-1:kfine)=y; zf(0:nf-1:kfine)=z
  do jk=1,kfine-1
     w=real(jk)/kfine; wc=1.-w
     xf(jk:jk+nmmkfine:kfine)=wc*x(0:nmm)+w*x(1:nm)
     yf(jk:jk+nmmkfine:kfine)=wc*y(0:nmm)+w*y(1:nm)
     zf(jk:jk+nmmkfine:kfine)=wc*z(0:nmm)+w*z(1:nm)
  enddo
end subroutine gglline3
!=============================================================================
!            SUBROUTINE GGFLINE1
! INTERPOLATE A SEGMENT OF N POINTS A TO FINER STRING AF
! USING SMOOTH PIECEWISE CUBIC METHOD WITH REFINEMENT FACTOR KFINE
!=============================================================================
subroutine ggfline1(n,nf,a,af)
  use pfstatus
  implicit none
  integer,intent(in):: n,nf
  real,dimension(0:n-1), intent(in) :: a
  real,dimension(0:nf-1),intent(out):: af
!-----------------------------------------------------------------------------
  integer:: jk,kf,nm1,nm2,nm3,nm4,nfm,nm3kf
  real   :: w0,w1,w2,w3
!=============================================================================
   if(kfine_a /= kfine)call ggsetw(kfine,fw1,fwn,fw); kfine_a=kfine
  nm1=n-1; nm2=n-2; nm3=n-3; nm4=n-4;nm3kf=nm3*kfine
  nfm=nm1*kfine; if(nfm /= nf-1)stop 'In ggfline; n, nf and kfine inconsistent'
  af(0:nfm:kfine)=a(0:nm1);
  do jk=1,kfine-1
     af(jk)    =dot_product(fw1(:,jk),a(0:2))
     af(nfm-jk)=dot_product(fw1(:,jk),a(nm1:nm3:-1))
     w0=fw(1,jk); w1=fw(2,jk); w2=fw(3,jk); w3=fw(4,jk)
     af(jk+kfine:jk+nm3kf:kfine)= &
          w0*a(0:nm4)+w1*a(1:nm3)+w2*a(2:nm2)+w3*a(3:nm1)
  enddo
end subroutine ggfline1
!=============================================================================
!            SUBROUTINE GGFLINE2
! INTERPOLATE A SEGMENT OF N POINTS X,Y TO FINER STRING XF,YF
! USING SMOOTH PIECEWISE CUBIC METHOD WITH REFINEMENT FACTOR KFINE
!=============================================================================
subroutine ggfline2(n,nf,x,y,xf,yf)
  use pfstatus
  implicit none
  integer,intent(in):: n,nf
  real,dimension(0:n-1), intent(in) :: x,y
  real,dimension(0:nf-1),intent(out):: xf,yf
!-----------------------------------------------------------------------------
  integer:: jk,kf,nm1,nm2,nm3,nm4,nfm,nm3kf
  real   :: w0,w1,w2,w3
!=============================================================================
   if(kfine_a /= kfine)call ggsetw(kfine,fw1,fwn,fw); kfine_a=kfine
  nm1=n-1; nm2=n-2; nm3=n-3; nm4=n-4;nm3kf=nm3*kfine
  nfm=nm1*kfine; if(nfm /= nf-1)stop 'In ggfline; n, nf and kfine inconsistent'
  xf(0:nfm:kfine)=x; yf(0:nfm:kfine)=y
  do jk=1,kfine-1
     xf(jk)    =dot_product(fw1(:,jk),x(0:2))
     yf(jk)    =dot_product(fw1(:,jk),y(0:2))
     xf(nfm-jk)=dot_product(fw1(:,jk),x(nm1:nm3:-1))
     yf(nfm-jk)=dot_product(fw1(:,jk),y(nm1:nm3:-1))
     w0=fw(1,jk); w1=fw(2,jk); w2=fw(3,jk); w3=fw(4,jk)
     xf(jk+kfine:jk+nm3kf:kfine)= &
          w0*x(0:nm4)+w1*x(1:nm3)+w2*x(2:nm2)+w3*x(3:nm1)
     yf(jk+kfine:jk+nm3kf:kfine)= &
          w0*y(0:nm4)+w1*y(1:nm3)+w2*y(2:nm2)+w3*y(3:nm1)
  enddo
end subroutine ggfline2

!=============================================================================
!            SUBROUTINE GGFLINE3
! INTERPOLATE A SEGMENT OF N POINTS X,Y,Z TO FINER STRING XF,YF,XF
! USING SMOOTH PIECEWISE CUBIC METHOD WITH REFINEMENT FACTOR KFINE
!=============================================================================
subroutine ggfline3(n,nf,x,y,z,xf,yf,zf)
  use pfstatus
  implicit none
  integer,intent(in):: n,nf
  real,dimension(0:n-1), intent(in) :: x,y,z
  real,dimension(0:nf-1),intent(out):: xf,yf,zf
!-----------------------------------------------------------------------------
  integer:: jk,kf,nm1,nm2,nm3,nm4,nfm,nm3kf
  real   :: w0,w1,w2,w3
!=============================================================================
   if(kfine_a /= kfine)call ggsetw(kfine,fw1,fwn,fw); kfine_a=kfine
  nm1=n-1; nm2=n-2; nm3=n-3; nm4=n-4;nm3kf=nm3*kfine
  nfm=nm1*kfine; if(nfm /= nf-1)stop 'In ggfline; n, nf and kfine inconsistent'
  xf(0:nfm:kfine)=x; yf(0:nfm:kfine)=y; zf(0:nfm:kfine)=z
  do jk=1,kfine-1
     xf(jk)    =dot_product(fw1(:,jk),x(0:2))
     yf(jk)    =dot_product(fw1(:,jk),y(0:2))
     zf(jk)    =dot_product(fw1(:,jk),z(0:2))
     xf(nfm-jk)=dot_product(fw1(:,jk),x(nm1:nm3:-1))
     yf(nfm-jk)=dot_product(fw1(:,jk),y(nm1:nm3:-1))
     zf(nfm-jk)=dot_product(fw1(:,jk),z(nm1:nm3:-1))
     w0=fw(1,jk); w1=fw(2,jk); w2=fw(3,jk); w3=fw(4,jk)
     xf(jk+kfine:jk+nm3kf:kfine)= &
          w0*x(0:nm4)+w1*x(1:nm3)+w2*x(2:nm2)+w3*x(3:nm1)
     yf(jk+kfine:jk+nm3kf:kfine)= &
          w0*y(0:nm4)+w1*y(1:nm3)+w2*y(2:nm2)+w3*y(3:nm1)
     zf(jk+kfine:jk+nm3kf:kfine)= &
          w0*z(0:nm4)+w1*z(1:nm3)+w2*z(2:nm2)+w3*z(3:nm1)
  enddo
end subroutine ggfline3

subroutine ggsetw(kx,xw1,xwn,xw)
  implicit none
  integer,intent(in):: kx
  real,intent(out):: xw1(2:4,4),xwn(3,4),xw(4,4)
  integer:: kxm,jkx
  real:: d,x,x1,xx
  kxm=kx-1
  d=1./kx
  do jkx=1,kxm
     x=d*jkx; x1=1.-x; xx=x*x
     xwn(1,jkx)=-x/2.+xx/2.;      xwn(2,jkx)=1.-xx;   xwn(3,jkx)=x/2.+xx/2.
     xw1(2,jkx)=1.-3.*x/2.+xx/2.; xw1(3,jkx)=2.*x-xx; xw1(4,jkx)=-x/2.+xx/2.
     xw(1,jkx)=x1*xwn(1,jkx)
     xw(2,jkx)=x1*xwn(2,jkx)+x*xw1(2,jkx)
     xw(3,jkx)=x1*xwn(3,jkx)+x*xw1(3,jkx)
     xw(4,jkx)=	             x*xw1(4,jkx)
  enddo
end subroutine ggsetw

subroutine gcaption(buf) ! write a centered caption
  use pfstatus
  implicit none
  character(len=*),intent(IN):: buf
  integer:: lfull,i
  real:: half
  if(.not.lopen)call ggopen
  call ggupdate; call ggupdatex
  lfull=0
  do i=1,80
     if(buf(i:i)=='%')exit
     lfull=i
  enddo
  half=lfull/2.
  half=wstring(buf,lfull)/2.
  call greltx(-half,0.)
  call groman(buf)
end subroutine gcaption

subroutine gxax(x1,x2,y)
  implicit none
  real,intent(in):: x1,x2,y
  real,parameter:: scaleh=99.096576,tscale=1.,ytn=.8,yt0=-2.5,ytp=-2.2,&
       ypc=-1.,ync=1.
  integer:: it1,j1,j2,k1,k2,idiv,kount,ix,i,j,k,itick,ilabel
  real:: rt1,div,div5,y_dev,x,y1,y2,y3,wchh,wchw,chh
  character*80 buft
  character*13 fedlst(6)
  data fedlst/'(f10.1,''%'')','(f10.2,''%'')','(f10.3,''%'')'&
             ,'(f10.4,''%'')','(f10.5,''%'')','(f10.6,''%'')'/
  call gginixax
  call gqtick(itick)
  call gqchh(chh)
  call gline(x1,y,x2,y) ! draw the axis
  call ggdecade(x1,x2,-1.,it1,rt1,div,div5,j1,j2,k1,k2)
  wchh=ggyttow(1.)
  y1=y+wchh*tscale*.5*(itick-1); 
  y2=y+wchh*tscale*.5*(itick+1)!  <- long ticks
  do j=j1,j2; x=j*div; call gline(x,y1,x,y2); enddo 
  y1=y+wchh*tscale*.25*(itick-1); 
  y2=y+wchh*tscale*.25*(itick+1)!  <- short ticks
  do k=k1,k2; if(mod(k,5)==0)cycle;x=k*div5;call gline(x,y1,x,y2); enddo
  call gqlabel(ilabel)
  y3=y
  if(ilabel==1)then
     select case(itick)
        case(-1); y3=y+ytn*wchh
        case(0);  y3=y+yt0*wchh
        case(1);  y3=y+ytp*wchh
     end select
     if(it1 >= 0)then 
        idiv=(div*1.00001)
        do j=j1,j2
           ix=idiv*j
           x=div*j
           write(unit=buft,fmt='(i7,''%'')')ix
           call ggsquash(buft,kount)
           call gswnxy(x,y3)
           call gcaption(buft)
        enddo
     else			   ! print numbers in f-format	   ! it1
        idiv=5
        if(rt1 <= 0.)idiv=1
        do j=j1,j2				 
           ix=idiv*j
           x=div*j
           if(ix == 0)then			 
              write(unit=buft,fmt='(1x,i1,''%'')')ix
           else					 
              if(it1 < -6)stop 'label out of bounds in gxax'
              write(unit=buft,fmt=fedlst(-it1))x
              call ggsquash(buft,kount)
           endif
           call gswnxy(x,y3)
           call gcaption(buft)
        enddo
     endif
  endif
  x=(x1+x2)*.5				
  if(itick >= 0)then
     y3=y3+ypc*wchh; call ggscapb(x,y3); else		  
     y3=y3+ync*wchh; call ggscapt(x,y3); 
  endif
  call gswnxy(x,y3)
  call ggfinxax
end subroutine gxax

subroutine gtxax(x1,x2,y,ggtform)
implicit none
real,intent(IN ):: x1,x2,y
real:: ggtform
!external:: ggtform
  real,parameter:: scaleh=99.096576,tscale=1.,ytn=.8,yt0=-2.5,ytp=-2.2,&
       ypc=-1.,ync=1.
  integer:: it1,j1,j2,k1,k2,idiv,kount,ix,i,j,k,itick,ilabel
  real:: rt1,div,div5,y_dev,x,y1,y2,y3,wchh,wchw,chh
  character*80 buft
  character*13 fedlst(6)
  data fedlst/'(f10.1,''%'')','(f10.2,''%'')','(f10.3,''%'')'&
             ,'(f10.4,''%'')','(f10.5,''%'')','(f10.6,''%'')'/
  call gginixax
  call gqtick(itick)
  call gqchh(chh)
  call gline(ggtform(x1),y,ggtform(x2),y) ! draw the axis
  call ggdecade(x1,x2,-1.,it1,rt1,div,div5,j1,j2,k1,k2)
  wchh=ggyttow(1.)
  y1=y+wchh*tscale*.5*(itick-1); 
  y2=y+wchh*tscale*.5*(itick+1)!  <- long ticks
  do j=j1,j2; x=j*div; call gline(ggtform(x),y1,ggtform(x),y2); enddo 
  y1=y+wchh*tscale*.25*(itick-1); 
  y2=y+wchh*tscale*.25*(itick+1)!  <- short ticks
  do k=k1,k2; if(mod(k,5)==0)cycle;x=k*div5;
     call gline(ggtform(x),y1,ggtform(x),y2); enddo
  call gqlabel(ilabel)
  y3=y
  if(ilabel==1)then
     select case(itick)
        case(-1); y3=y+ytn*wchh
        case(0);  y3=y+yt0*wchh
        case(1);  y3=y+ytp*wchh
     end select
     if(it1 >= 0)then 
        idiv=(div*1.00001)
        do j=j1,j2
           ix=idiv*j
           x=div*j
           write(unit=buft,fmt='(i7,''%'')')ix
           call ggsquash(buft,kount)
           call gswnxy(ggtform(x),y3)
           call gcaption(buft)
        enddo
     else			   ! print numbers in f-format	   ! it1
        idiv=5
        if(rt1 <= 0.)idiv=1
        do j=j1,j2				 
           ix=idiv*j
           x=div*j
           if(ix == 0)then			 
              write(unit=buft,fmt='(1x,i1,''%'')')ix
           else					 
              if(it1 < -6)stop'label out of bounds in gtxax'
              write(unit=buft,fmt=fedlst(-it1))x
              call ggsquash(buft,kount)
           endif
           call gswnxy(ggtform(x),y3)
           call gcaption(buft)
        enddo
     endif
  endif
  x=(ggtform(x1)+ggtform(x2))*.5				
  if(itick >= 0)then
     y3=y3+ypc*wchh; call ggscapb(x,y3); else		  
     y3=y3+ync*wchh; call ggscapt(x,y3); 
  endif
  call gswnxy(x,y3)
  call ggfinxax
end subroutine gtxax

subroutine gyax(y1,y2,x)
  implicit none
  real,parameter:: scaleh=99.096576,tscale=1.,xtn=.7,xt0=.0,xtp=.2,&
       xpc=-1.5,xnc=1.
  real,intent(in):: y1,y2,x
  integer:: itick,mcount,it1,j1,j2,k1,k2,j,iy,kount,k,ilabel,idiv
  real:: chh,rt1,div,div5,wchh,wchw,x1,x2,x3,y,wcount,wcountm
  character*80 buft
  character*13 fedlst(6)
  data fedlst/'(f10.1,''%'')','(f10.2,''%'')','(f10.3,''%'')' &
       ,'(f10.4,''%'')','(f10.5,''%'')','(f10.6,''%'')'/
  call gginixax
  call gqtick(itick)
  call gqchh(chh)
  call gline(x,y1,x,y2)
  call ggdecade(y1,y2,-1.,it1,rt1,div,div5,j1,j2,k1,k2)
  wchh=ggyttow(1.); wchw=ggxttow(1.)
  x1=x+wchw*tscale*.5*(itick-1); 
  x2=x+wchw*tscale*.5*(itick+1)!  <- long ticks
  do j=j1,j2; y=j*div; call gline(x1,y,x2,y); enddo 
  x1=x+wchw*tscale*.25*(itick-1); 
  x2=x+wchw*tscale*.25*(itick+1)!  <- short ticks
  do k=k1,k2; if(mod(k,5)==0)cycle;y=k*div5;call gline(x1,y,x2,y); enddo
  call gqlabel(ilabel)
  wcountm=0.
  x1=x
  if(ilabel==1)then
     select case(itick)
        case(-1); x1=x+xtn*wchw
        case(0);  x1=x+xt0*wchw
        case(1);  x1=x+xtp*wchw
     end select
     if(it1 >= 0)then 
        idiv=(div*1.00001)
        do j=j1,j2
           iy=idiv*j
           y=div*j
           write(unit=buft,fmt='(i7,''%'')')iy
           call ggsquash(buft,kount)
           wcount=wstring(buft,kount)
           wcountm=max(wcount,wcountm)
           x3=x1
           if(itick >= 0)x3=x1-(wcount+1.)*wchw
           y=y-.5*wchh 
           call gswnxy(x3,y)
           call groman(buft)
        enddo
     else			  ! numbers in f-format
        idiv=5; if(rt1 <= 0.)idiv=1
        do j=j1,j2
           y=div*j
           iy=idiv*j
           if(iy==0)then
              write(unit=buft,fmt='(i1,''%'')')iy
              wcount=1.
            else
              if(it1 < -6)stop 'label out of bounds in gyax'
              write(unit=buft,fmt=fedlst(-it1))y
              call ggsquash(buft,kount)
              wcount=wstring(buft,kount)
              wcountm=max(wcountm,wcount)
           endif
           x3=x1
           if(itick >= 0)x3=x1-(wcount+1.)*wchw
           y=y-.5*wchh 
           call gswnxy(x3,y)
           call groman(buft)
        enddo
     endif
  endif
  y=(y1+y2)*.5
  if(itick >= 0)then; 
     x3=x1+(xpc-wcountm)*wchw; call ggscapl(x3,y); else
     x3=x1+(xnc+wcountm)*wchw; call ggscapr(x3,y); 
  endif
  call gswnxy(x3,y)
  call ggfinxax
end subroutine gyax

subroutine gtyax(y1,y2,x,ggtform)
  implicit none
  real,parameter:: scaleh=99.096576,tscale=1.,xtn=.7,xt0=.0,xtp=.2,&
       xpc=-1.5,xnc=1.
  real,intent(in):: y1,y2,x
!  external:: ggtform
  real:: ggtform
  integer:: itick,mcount,it1,j1,j2,k1,k2,j,iy,kount,k,ilabel,idiv
  real:: chh,rt1,div,div5,wchh,wchw,x1,x2,x3,y,wcount,wcountm
  character*80 buft
  character*13 fedlst(6)
  data fedlst/'(f10.1,''%'')','(f10.2,''%'')','(f10.3,''%'')' &
       ,'(f10.4,''%'')','(f10.5,''%'')','(f10.6,''%'')'/
  call gginixax
  call gqtick(itick)
  call gqchh(chh)
  call gline(x,ggtform(y1),x,ggtform(y2))
  call ggdecade(y1,y2,-1.,it1,rt1,div,div5,j1,j2,k1,k2)
  wchh=ggyttow(1.); wchw=ggxttow(1.)
  x1=x+wchw*tscale*.5*(itick-1); 
  x2=x+wchw*tscale*.5*(itick+1)!  <- long ticks
  do j=j1,j2; y=j*div; call gline(x1,ggtform(y),x2,ggtform(y)); enddo 
  x1=x+wchw*tscale*.25*(itick-1); 
  x2=x+wchw*tscale*.25*(itick+1)!  <- short ticks
  do k=k1,k2; if(mod(k,5)==0)cycle;y=k*div5;
     call gline(x1,ggtform(y),x2,ggtform(y)); enddo
  call gqlabel(ilabel)
  wcountm=0.
  x1=x
  if(ilabel==1)then
     select case(itick)
        case(-1); x1=x+xtn*wchw
        case(0);  x1=x+xt0*wchw
        case(1);  x1=x+xtp*wchw
     end select
     if(it1 >= 0)then 
        idiv=(div*1.00001)
        do j=j1,j2
           iy=idiv*j
           y=div*j
           write(unit=buft,fmt='(i7,''%'')')iy
           call ggsquash(buft,kount)
           wcount=wstring(buft,kount)
           wcountm=max(wcount,wcountm)
           x3=x1
           if(itick >= 0)x3=x1-(wcount+1.)*wchw
           y=y-.5*wchh 
           call gswnxy(x3,ggtform(y))
           call groman(buft)
        enddo
     else			  ! numbers in f-format
        idiv=5; if(rt1 <= 0.)idiv=1
        do j=j1,j2
           y=div*j
           iy=idiv*j
           if(iy==0)then
              write(unit=buft,fmt='(i1,''%'')')iy
              wcount=1.
            else
              if(it1 < -6)stop 'label out of bounds in gyax'
              write(unit=buft,fmt=fedlst(-it1))y
              call ggsquash(buft,kount)
              wcount=wstring(buft,kount)
              wcountm=max(wcountm,wcount)
           endif
           x3=x1
           if(itick >= 0)x3=x1-(wcount+1.)*wchw
           y=y-.5*wchh 
           call gswnxy(x3,ggtform(y))
           call groman(buft)
        enddo
     endif
  endif
  y=(ggtform(y1)+ggtform(y2))*.5
  if(itick >= 0)then; 
     x3=x1+(xpc-wcountm)*wchw; call ggscapl(x3,y); else
     x3=x1+(xnc+wcountm)*wchw; call ggscapr(x3,y); 
  endif
  call gswnxy(x3,y)
  call ggfinxax
end subroutine gtyax

subroutine gxlax(x1s,x2s,y) ! plot logarithmic x axis ######################
  implicit none
  real,intent(in):: x1s,x2s,y
  real,parameter:: tscale=1.,ypc=-2.5,ync=2.1
  character*80 buft
  integer:: itick,j1,j2,j,k,k1,k2,i,ilabel,kount
  real:: x1,x2,chh,dchh,dchhh,dchhq,d,dp,xt1,xt2,ycc,y1,y2,x,y1c,y2c,&
       rk,alk,x3c,y3c,x3,y3,xcapb,ycapb,xcapt,ycapt,wchh
  x1=x1s
  x2=x2s
  y3=y
  call gginixax
  call gqtick(itick)
  call gqchh(chh)
  dchh=99.096576*chh
  dchhh=dchh*.5
  dchhq=dchh*.25
  call gline(x1s,y,x2s,y)
  wchh=ggyttow(1.)
  d=x2-x1
  dp=d*1.00001
  xt1=x2-dp
  xt2=x1+dp
  j1=nint(xt1+.5)
  j2=nint(xt2-.5)
  y3c=ggywtod(y)
  ycc=y3c+dchhh*itick
  y1=ggydtow(ycc-dchhh)
  y2=ggydtow(ycc+dchhh)
  y1=y+wchh*tscale*.5*(itick-1)
  y2=y+wchh*tscale*.5*(itick+1)

!  long ticks:
  do j=j1,j2; x=j; call gline(x,y1,x,y2); enddo
  ycc=y3c+dchhq*itick
  y1c=ycc-dchhq
  y2c=ycc+dchhq
  y1=ggydtow(y1c)
  y2=ggydtow(y2c)
  y1=y+wchh*tscale*.25*(itick-1); 
  y2=y+wchh*tscale*.25*(itick+1)!  <- short ticks

!  short ticks:
  do k=2,9              
     rk=k
     alk=alog10(rk)
     k1=nint(x1-alk+.5)
     k2=nint(x2-alk-.5)
     do i=k1,k2; x=alk+i; call gline(x,y1,x,y2); enddo
  enddo
  call gqlabel(ilabel)
  if(ilabel==1)then
     if(itick==1)then
        y3c=y1c-.8*dchh
     elseif(itick==0)then
        y3c=y1c-1.*dchh
     elseif(itick==-1)then
        y3c=y2c+1.8*dchh
     endif
  else
     y3c=y3c-.3*dchh
  endif
  if(ilabel==1)then
     do j=j1,j2
        x=j
        x3c=ggxwtod(x)
        call gsdvxy(x3c,y3c)
        call greltx(-1.3,-1.4)
        call groman('10%')
        call gsup
        write(unit=buft,fmt='(i7,''%'')')j
        call ggsquash(buft,kount)
        call groman(buft)
        call gunsup
     enddo
  endif
  x=(x1s+x2s)*.5	
  if(itick >= 0)then
     y3=y3+ypc*wchh; call ggscapb(x,y3); else
     y3=y3+ync*wchh; call ggscapt(x,y3);
  endif
  call gswnxy(x,y3)
  call ggfinxax
end subroutine gxlax

subroutine gylax(y1s,y2s,x) ! plot logarithmic y-axis ########################
  real,intent(in):: y1s,y2s,x
  character*80 buft
  y1=y1s
  y2=y2s
  call gginixax
  call gqtick(itick)
  call gqchh(chh)
  dchh=99.096576*chh
  dchhh=dchh*.5
  dchhq=dchh*.25
  dchw=dchhh*150./81.28
  call gline(x,y1s,x,y2s)
  d=y2-y1
  dp=d*1.00001
  yt1=y2-dp
  yt2=y1+dp
  j1=nint(yt1+.5)
  j2=nint(yt2-.5)
  x3c=ggxwtod(x)
  xcc=x3c+dchhh*itick
  x1=ggxdtow(xcc-dchhh)
  x2=ggxdtow(xcc+dchhh)
!  long ticks:
  do j=j1,j2            
     y=j
     call gline(x1,y,x2,y)
  enddo
  xcc=x3c+dchhq*itick
  x1c=xcc-dchhq
  x2c=xcc+dchhq
  x1=ggxdtow(x1c)
  x2=ggxdtow(x2c)
!  short ticks:
  do k=2,9
     rk=k
     alk=alog10(rk)
     k1=nint(y1-alk+.5)
     k2=nint(y2-alk-.5)
     do i=k1,k2
        y=alk+i
        call gline(x1,y,x2,y)
     enddo
  enddo
  if(itick >= 0)then
     x3c=x1c-2.*dchw
  else
     x3c=x1c+3.*dchw
  endif
  call gqlabel(ilabel)
  if(ilabel == 1)then
     do j=j1,j2
        y=j
        y3c=ggywtod(y)
        call gsdvxy(x3c,y3c)
        call greltx(-2.,-.5)
        call groman('10%')
        call gsup
        write(unit=buft,fmt='(i7,''%'')')j
        call ggsquash(buft,kount)
        call groman(buft)
        call gunsup
     enddo
  endif
  y3=(y1+y2)*.5						 
  y3c=ggywtod(y3)
  if(itick >= 0)then
     x3c=x3c-2.2*dchw
     x3=ggxdtow(x3c)
     xcapl=x3
     ycapl=y3
     call ggscapl(xcapl,ycapl)
  else
     x3c=x3c+1.7*dchw
     x3=ggxdtow(x3c)
     xcapr=x3
     ycapr=y3
     call ggscapr(xcapr,ycapr)
  endif
  call gsdvxy(x3c,y3c)
  call ggfinxax
end subroutine gylax

!============================================================================
subroutine garr(u,v,dt2)
implicit none
real,dimension(:,:),  intent(IN):: u,v
real,                 intent(IN):: dt2
!----------------------------------------------------------------------------
real                :: x0,x1,y0,y1,dt,u0,v0,ux,uy,ut,vx,vy,vt,ex,ey
real,dimension(-1:1):: x,y,z
integer             :: nx,ny,i,j,fsave
!============================================================================
nx=size(u,1); ny=size(u,2)
if(nx/=size(v,1) .or. ny/=size(v,2))stop 'In garr; incompatible dimensions'
dt=dt2/2
call gqfine(fsave)
call gsfine(5)
call gqnt(j); 
call gqwn(j,x0,x1,y0,y1)
call grect(x0,x1,y0,y1)
call gspli(0)
ex=(x1-x0)/(nx-1)
ey=(y1-y0)/(ny-1)
do j=1,ny
do i=1,nx
   x(0)=x0+(i-1)*ex
   y(0)=y0+(j-1)*ey
   z=1.
   u0=u(i,j)
   v0=v(i,j)
   if(i==1)then
      ux=(u(i+1,j)-u(i,j))/ex
      vx=(v(i+1,j)-v(i,j))/ex
   elseif(i==nx)then
      ux=(u(i,j)-u(i-1,j))/ex
      vx=(v(i,j)-v(i-1,j))/ex
   else
      ux=(u(i+1,j)-u(i-1,j))/(2*ex)
      vx=(v(i+1,j)-v(i-1,j))/(2*ex)
   endif
   if(j==1)then
      uy=(u(i,j+1)-u(i,j))/ey
      vy=(v(i,j+1)-v(i,j))/ey
   elseif(j==ny)then
      uy=(u(i,j)-u(i,j-1))/ey
      vy=(v(i,j)-v(i,j-1))/ey
   else
      uy=(u(i,j+1)-u(i,j-1))/(2*ey)
      vy=(v(i,j+1)-v(i,j-1))/(2*ey)
   endif
   ut=ux*u0+uy*v0
   vt=vx*u0+vy*v0
   x(1) =x(0)+(u0+ut*dt/2)*dt
   y(1) =y(0)+(v0+vt*dt/2)*dt
   x(-1)=x(0)-(u0-ut*dt/2)*dt
   y(-1)=y(0)-(v0-vt*dt/2)*dt
   if(x(1)<x0 .or. x(1)>x1 .or. y(1)<y0 .or. y(1)>y1)z(1)=-1
   if(x(-1)<x0 .or. x(-1)>x1 .or. y(-1)<y0 .or. y(-1)>y1)z(-1)=-1
   call gpl(3,x,y,z)
enddo
enddo
call gsfine(1)
dt=dt/5
do j=1,ny
do i=1,nx
   x(0)=x0+(i-1)*ex
   y(0)=y0+(j-1)*ey
   u0=u(i,j)
   v0=v(i,j)
   x(-1)=x(0)-u0*dt+v0*dt
   y(-1)=y(0)-v0*dt-u0*dt
   x( 1)=x(0)-u0*dt-v0*dt
   y( 1)=y(0)-v0*dt+u0*dt
   x( 0)=x(0)+u0*dt
   y( 0)=y(0)+v0*dt
   call gpl(3,x,y)
enddo
enddo
call gsfine(fsave)
end subroutine garr

subroutine gdcade(ad,mcon,dinck,mode)
  implicit none
  real,   intent(IN) :: ad
  integer,intent(IN) :: mcon,mode
  real,   intent(OUT):: dinck
  integer:: i1
  real:: a2,a4,a5,bd,r1,t1,acon,b1
  data a2/.30103/,a4/.60206/,a5/.69897/
  acon=real(mcon)-1.
  bd=alog10(ad/acon)
  if(mode==2)then
     bd=bd+a5
  elseif(mode==0)then
     bd=bd+a2
  endif
  i1=nint(bd-.5)
  b1=i1
  r1=bd-i1
  t1=10**b1
  if(mode==2)then
     if(r1<=a5)then
        dinck=t1
     else
        dinck=2*t1
     endif
  else
     if(r1<=a2)then
        dinck=t1
     elseif(r1<=a4)then
        dinck=2*t1
     else
        dinck=5*t1
     endif
  endif
end subroutine gdcade

subroutine gstang(xxrun,yyrise)
  use pfstatus
  if(.not.lopen)call ggopen
  s=sqrt(xxrun**2+yyrise**2)
  xrun=xxrun/s
  yrise=yyrise/s
  return
  entry gqtang(xxrun,yyrise)
  xxrun=xrun
  yyrise=yrise
end subroutine gstang

subroutine gsbdiv(ibdiv)
  use pfstatus
  nbdiv=ibdiv
  nbdivp=nbdiv+1
  return
entry gqbdiv(ibdiv)
ibdiv=nbdiv
end subroutine gsbdiv

subroutine gbordr(x1,y1,x2,y2,ggtform)
  use pfstatus
  external ggtform
  if(jv /= 0)call ggenlin
  call gginixax
  difx=(x2-x1)/nbdiv; dify=(y2-y1)/nbdiv
  do i=1,nbdivp; xarr(i)=x1+(i-1)*difx; yarr(i)=y1+(i-1)*dify; enddo
  call gpl(nbdivp,xarr,yarr,ggtform)
  call ggfinxax
end subroutine gbordr

subroutine gssty(isty)
  use pfstatus
                     ksty=isty;   return
entry gssdar(asdar); rsdar=asdar; return
entry gsspen(jspen); ispen=jspen; return
entry gssodd(asodd); rsodd=asodd; return
entry gdsty;         ksty=1;      return
entry gdsdar;        rsdar=.5;    return
entry gdspen;        ispen=1;     return
entry gdsodd;        rsodd=0.;    return
entry gqsty(isty);   isty=ksty;   return
entry gqsdar(asdar); asdar=rsdar; return
entry gqspen(jspen); jspen=ispen; return
entry gqsodd(asodd); asodd=rsodd
end subroutine gssty

subroutine gxshade(nb,x,y,ysep)
! apply cross-shading to the polygonal loop (x,y) by invoking horizontal
! and vertical shading (without staggering the dashes for broken line styles)
  dimension x(nb),y(nb)
  call gdsodd
  call ghshade(nb,x,y,ysep)
  call gvshade(nb,x,y,ysep)
end subroutine gxshade

subroutine gwshade(nb,x,y,ysep)
! apply a "weave"-style shading by invoking broken horizontal and vertical
! shading with dashes staggered through a setting of rsodd=.5 on the vertical
  dimension x(nb),y(nb)
  call gqsty(kstys)
  call gssty(2)
  call gdsodd
  call ghshade(nb,x,y,ysep)
  call gssodd(.5)
  call gvshade(nb,x,y,ysep)
  call gdsodd
  call gssty(kstys)
end subroutine gwshade

subroutine ghshade(nb,x,y,ysep)
!  apply horizontal shading to the interior of the polygonal loop, (x,y)
  use pfstatus
  dimension x(nb),y(nb)
  data sep0/64./,rh/.7071068/
entry gshade(nb,x,y,ysep)
  rxx=bx/(ysep*sep0)
  ryx=0.
  rxy=0.
  ryy=by/(ysep*sep0)
  call ggshade(nb,x,y,rxx,ryx,rxy,ryy)
  return

entry gvshade(nb,x,y,ysep)
!  apply vertical shading to the interior of the polygonal loop, (x,y)
  ryx=bx/(ysep*sep0)
  rxy=-by/(ysep*sep0)
  rxx=0.
  ryy=0.
  call ggshade(nb,x,y,rxx,ryx,rxy,ryy)
  return

entry grshade(nb,x,y,ysep)
!  apply right-slanted shading to the interior of the polygonal loop, (x,y)
  rxx=bx/(ysep*sep0)
  ryy=by/(ysep*sep0)
  rxy=rh*ryy
  ryy=rh*ryy
  ryx=-rh*rxx
  rxx=rh*rxx
  call ggshade(nb,x,y,rxx,ryx,rxy,ryy)
  return

entry glshade(nb,x,y,ysep)
!  apply left-slanted shading to the interior of the polygonal loop, (x,y)
  rxx=bx/(ysep*sep0)
  ryy=by/(ysep*sep0)
  rxy=-rh*ryy
  ryy=rh*ryy
  ryx=rh*rxx
  rxx=rh*rxx
  call ggshade(nb,x,y,rxx,ryx,rxy,ryy)
  return
end subroutine ghshade

subroutine ggshade(nb,x,y,rxx,ryx,rxy,ryy)
  use pfstatus
  integer,           intent(IN):: nb
  real,dimension(nb),intent(IN):: x,y
  real,              intent(IN):: rxx,ryx,rxy,ryy
  real,dimension(20):: xcut
  logical la,lb
  data eps/.0001/
  if(.not.lopen)call ggopen
  call ggupdate
  if(ilpen /= ispen)then
     ilpen=ispen
     wid=pwid(ilpen)
     wid_a=wid
     call ggfmt107(ggunit,wid)
  endif

  d=1./(rxx*ryy-ryx*rxy)
  qxx=bx*ryy*d
  qyx=-by*ryx*d
  qxy=-bx*rxy*d
  qyy=by*rxx*d
  if(nb <= 2)return
  ymin=ryx*x(1)+ryy*y(1)
  ymax=ymin
  do ib=2,nb
     yb=ryx*x(ib)+ryy*y(ib)
     ymin=min(ymin,yb)
     ymax=max(ymax,yb)
  enddo
  do iline=nint(ymin+.5),nint(ymax-.5)
     iline2=mod(iline,2)
     yh=iline
     ncut=0
     ia=nb
     xia=rxx*x(ia)+rxy*y(ia)
     yia=ryx*x(ia)+ryy*y(ia)
     la=yia.ge.yh
     do ib=1,nb
	xib=rxx*x(ib)+rxy*y(ib)
	yib=ryx*x(ib)+ryy*y(ib)
	lb=yib.ge.yh
        if((la .and. .not.lb) .or. (lb .and. .not.la))then
           if(ncut==20)stop 'region too complex to allow shading'
           ncut=ncut+1
           ya=yia-yh
           yb=yib-yh
           xcut(ncut)=(xia*yb-xib*ya)/(yb-ya)
	endif
	la=lb
	ia=ib
	xia=xib
	yia=yib
     enddo
     do icut=1,ncut-1
	xkcut=xcut(icut)
	kcut=icut
	do jcut=icut+1,ncut
           if(xcut(jcut).lt.xkcut)then
              kcut=jcut
              xkcut=xcut(jcut)
           endif
	enddo
	if(kcut /= icut)then
           xcut(kcut)=xcut(icut)
           xcut(icut)=xkcut
	endif
     enddo
     icut1=1
400  icut2=icut1+1
     if(icut2 <= ncut)then
	if(xcut(icut2)-xcut(icut1) < eps)then
           do icut=icut2+1,ncut
              xcut(icut-2)=xcut(icut)
           enddo
           ncut=ncut-2
	else
           icut1=icut2
	endif
	goto 400
     endif
     do ib=2,ncut,2
        ia=ib-1
        xa=xcut(ia)
        xb=xcut(ib)
        xc1=ax+qxx*xa+qxy*yh
        xc2=ax+qxx*xb+qxy*yh
        yc1=ay+qyx*xa+qyy*yh
        yc2=ay+qyx*xb+qyy*yh
        
        if(ksty==1)then
           call ggfmt106(ggunit,xc1,yc1,xc2,yc2)
           call ggfmt108(ggunit)
        elseif(ksty==2)then
           odd=iline2/2.+rsodd-rsdar/2.
           ic1=nint(xa/2.+.5-odd)
           ic2=nint(xb/2.-.5-odd)
           xc=((ic1-1)+odd)*2.
           xd=((ic1-1)+odd+rsdar)*2.
           xc=amax1(xc,xa)
           xd=amin1(xd,xb)
           if(xd > xc)then ! leading partial dash
              xc3=ax+qxx*xc+qxy*yh
              yc3=ay+qyx*xc+qyy*yh
              xc4=ax+qxx*xd+qxy*yh
              yc4=ay+qyx*xd+qyy*yh
              call ggfmt106(ggunit,xc3,yc3,xc4,yc4)
              call ggfmt108(ggunit)
           endif

           do ic=ic1,ic2-1	! all the complete dashes in the middle
              xc=(ic+odd)*2.
              xd=(ic+odd+rsdar)*2.
              xc3=ax+qxx*xc+qxy*yh
              yc3=ay+qyx*xc+qyy*yh
              xc4=ax+qxx*xd+qxy*yh
              yc4=ay+qyx*xd+qyy*yh
              call ggfmt106(ggunit,xc3,yc3,xc4,yc4)
              call ggfmt108(ggunit)
           enddo

           xc=((ic2)+odd)*2.
           xd=((ic2)+odd+rsdar)*2.
           xc=amax1(xc,xa)
           xd=amin1(xd,xb)
           if(xd.gt.xc)then ! trailing partial dash
              xc3=ax+qxx*xc+qxy*yh
              yc3=ay+qyx*xc+qyy*yh
              xc4=ax+qxx*xd+qxy*yh
              yc4=ay+qyx*xd+qyy*yh
              call ggfmt106(ggunit,xc3,yc3,xc4,yc4)
              call ggfmt108(ggunit)
           endif
        endif
     enddo
  enddo
end subroutine ggshade

subroutine gspli(jltp)	 ! set polyline index
  use pfconsts
  use pfstatus
  if(jltp<0 .or. jltp>nltp)stop 'In gspli; jltp out of bounds'
  if(lstroke)call ggfmt108(ggunit)
  iltp=jltp; ilpen=lpeni(iltp); iltype=ltypei(iltp); return
entry gqpli(jltp); jltp=iltp; return
entry gspmi(jmtp)
if(jmtp<0 .or. jmtp>nmtp)stop 'In gspmi; jmtp out of bounds'
imtp=jmtp; imtype=mtypei(imtp); impen=mpeni(imtp);
rpm=rpmi*45.; rc30=rpm*c30; rs30=rpm*s30
end subroutine gspli

subroutine ggreek(buf)
  use pfconsts
  use pfstatus
  implicit none
  character(len=*):: buf
  character*80 :: buf2
  character*1 ch
  integer:: ilpen_save,nich,nblank,ich,ifont,i,kchar
  real:: xc1,yc1,xdev
  if(.not.lopen)call ggopen
  ilpen_save=ilpen
  ilpen=itpen
  call ggupdate
  call gqfont(ifont)
  call gsfont(7)
  do ich=1,80
     ch=buf(ich:ich)
     kchar=ichar(ch)
     if(kchar<32.or.kchar>=127)stop 'character not recognized in subr. ggreek'
     if(ch=='%')goto 375
     nich=ich
  enddo
375 nblank=0
  do ich=nich,1,-1
     if(buf(ich:ich)/=' ')goto 376
     nblank=nblank+1
  enddo
376 if(nblank/=nich)then
     write(unit=buf2,fmt='(80a1)')(buf(i:i),i=1,nich-nblank)
     call gqwnxy(xc1,yc1)
     call gtx(xc1,yc1,buf2)
  endif
  xdev=nblank*wblank()
  call greltx(xdev,0.)
!  ifonts=ifont  ! return current font index to its former value
  call gsfont(ifont)
!  call ggfmt112(chh*cfudge)
  ilpen=ilpen_save
  contains
    real function wblank()
      character*80 buf
      buf(1:1)=' '
      wblank=wstring(buf,1)
    end function wblank
end subroutine ggreek

subroutine groman(buf)
  use pfconsts
  use pfstatus
  character(len=*):: buf
  character*80 buf2
  character*1 ch
  integer:: ilpen_save,kchar
  if(.not.lopen)call ggopen
  ilpen_save=ilpen
  ilpen=itpen
  call ggupdate
  do ich=1,80
     ch=buf(ich:ich)
     kchar=ichar(ch)
     if(kchar<32.or.kchar>=127)stop 'character not recognized in subr. groman'
     if(ch=='%')goto 375
     nich=ich
  enddo
375 nblank=0
  do ich=nich,1,-1
     if(buf(ich:ich)/=' ')goto 376
     nblank=nblank+1
  enddo
376 if(nblank/=nich)then
     write(unit=buf2,fmt='(80a1)')(buf(i:i),i=1,nich-nblank)
     call gqwnxy(xc1,yc1)
     call gtx(xc1,yc1,buf2)
  endif
  xdev=nblank*wblank()
  call greltx(xdev,0.)
  ilpen=ilpen_save
  contains
    real function wblank()
      character*80 buf
      buf(1:1)=' '
      wblank=wstring(buf,1)
    end function wblank
end subroutine groman

subroutine gswnxy(x,y)
  use pfstatus
                  if(.not.lopen)call ggopen; xc2=ax+bx*x; yc2=ay+by*y; return
entry gqwnxy(x,y);if(.not.lopen)call ggopen;x=(xc2-ax)/bx;y=(yc2-ay)/by;return
entry gsvpxy(x,y);if(.not.lopen)call ggopen; xc2=x*dy; yc2=y*dy;
call ggfmt100(ggunit,xc2,yc2); return
entry gqvpxy(x,y);if(.not.lopen)call ggopen; x=xc2/dy; y=yc2/dy; return
entry gsdvxy(x,y);if(.not.lopen)call ggopen; xc2=x; yc2=y
call ggfmt100(ggunit,xc2,yc2); return
entry gqdvxy(x,y);if(.not.lopen)call ggopen; x=xc2; y=yc2; return
end subroutine gswnxy

subroutine ghat	 ! put a caret over previous character
data up/.5/
call greltx(-1.,up); 
call groman('^%'); call greltx(0.,-up); return
entry gdash	 ! put a prime against previous character
call gsup; call gsub; 
call groman('/%'); call gunsub; call gunsup
end subroutine ghat

subroutine gsup
  use pfstatus
  data crat/.6/	! ratio of sizes of superscripts to standard chars.
  data txstep/.15/,tysup/1.3/,tysub/-.5/
                   call greltx(-txstep,tysup); chh=chh*crat
  call gschh(chh); call greltx(txstep,0.);                  return
entry gunsup;      call greltx(-txstep,0.); chh=chh/crat
  call gschh(chh); call greltx(txstep,-tysup);              return
entry gsub;        call greltx(-txstep,tysub); chh=chh*crat
  call gschh(chh); call greltx(txstep,0.);                  return
entry gunsub;      call greltx(-txstep,0.); chh=chh/crat
 call gschh(chh);  call greltx(txstep,-tysub);              return
end subroutine gsup

subroutine greltx(x,y)	! move pen in relative text coordinates
  use pfconsts
  use pfstatus
  if(.not.lopen)call ggopen
  xc2=xc2+(x*dchw*xrun-y*dchh*yrise)
  yc2=yc2+(x*dchw*yrise+y*dchh*xrun)
end subroutine greltx

subroutine gdcons
  use pfstatus
  ncon=0
end subroutine gdcons

subroutine gsctyp(n,jtyp) ! set contour types
  use pfconsts
  use pfstatus
  integer             :: n
  integer,dimension(n):: jtyp
  if(n>ncon_max)stop 'error in gsctyp. ncon > ncon_max'
  if(n<0)  stop 'error in gsctyp. ncon < 0'
                        ncon=n; ctyp(1:n)=jtyp(1:n); return
  entry gqctyp(n,jtyp); n=ncon; jtyp(1:n)=ctyp(1:n)
end subroutine gsctyp

subroutine gscval(n,aval) ! set contour values, types
  use pfconsts
  use pfstatus
  implicit none
  integer          :: n
  real,dimension(n):: aval
  if(n>ncon_max)stop 'error in gscval. ncon > ncon_max'
  if(n<0)  stop 'error in gscval. ncon < 0'
  lcval=.true.;         ncon=n; cval(1:n)=aval(1:n); return
entry gqcval(n,aval); n=0 
  if(lcval)then; n=ncon; aval(1:n)=cval(1:n); endif; return
entry gdcval; lcval=.false.; ncon=0; return
end subroutine gscval
  

subroutine ggscapxy(xcb,ycb,xct,yct,xcl,ycl,xcr,ycr)
call ggscapb(xcb,ycb); call ggscapt(xct,yct)
call ggscapl(xcl,ycl); call ggscapr(xcr,ycr)
end subroutine ggscapxy

subroutine gguform(x1,y1,x2,y2,z2);
implicit none;real:: x1,y1,x2,y2,z2; x2=x1;y2=y1;z2=1.;end subroutine gguform

subroutine ggfmt98(ggunit)
  integer,intent(IN):: ggunit
  write(ggunit,98)
98 FORMAT('%!PS'/                                                   &
        '/hpu {0.0711 mul } bind def'/'/hpc {10.66 mul } bind def'/ &
        ' % ggfmt98')
end subroutine ggfmt98

subroutine ggfmt99(ggunit)
  integer,intent(IN):: ggunit
  write(ggunit,99)
99 FORMAT('%!PS'/                                                   &
        '/orix {580} bind def /oriy {30} bind def'/                 &
        'orix oriy translate 90 rotate'/                            &
        '/hpu {0.0711 mul } bind def'/'/hpc {10.66 mul } bind def'/ &
        ' % ggfmt99')
end subroutine ggfmt99

subroutine ggfmt100(ggunit,x1,y1)
  integer,intent(IN):: ggunit
  real,intent(IN):: x1,y1
  write(ggunit,100)x1,y1
100 format(2(f7.1,' hpu '),'moveto % ggfmt100')
end subroutine ggfmt100

subroutine ggfmt101(ggunit,x1,y1)
  integer,intent(IN):: ggunit
  real,   intent(IN):: x1,y1
  write(ggunit,101)x1,y1
101 format(2(f7.1,' hpu '),'lineto   % ggfmt101')
end subroutine ggfmt101

subroutine ggfmt102(ggunit,x1,y1,x2,y2)
  integer,intent(IN):: ggunit
  real,intent(IN):: x1,y1,x2,y2
  call ggfmt101(ggunit,x1,y1)
  call ggfmt101(ggunit,x2,y2)
end subroutine ggfmt102

subroutine ggfmt103(ggunit,x1,y1,x2,y2,x3,y3)
  integer,intent(IN):: ggunit
  real,intent(IN):: x1,y1,x2,y2,x3,y3
  call ggfmt101(ggunit,x1,y1)
  call ggfmt101(ggunit,x2,y2)
  call ggfmt101(ggunit,x3,y3)
end subroutine ggfmt103

subroutine ggfmt104(ggunit,x1,y1,x2,y2,x3,y3,x4,y4)
  integer,intent(IN):: ggunit
  real,   intent(IN):: x1,y1,x2,y2,x3,y3,x4,y4
  call ggfmt101(ggunit,x1,y1)
  call ggfmt101(ggunit,x2,y2)
  call ggfmt101(ggunit,x3,y3)
  call ggfmt101(ggunit,x4,y4)
end subroutine ggfmt104

subroutine ggfmt105(ggunit,x1,y1,x2,y2,x3,y3,x4,y4)
  integer,intent(IN):: ggunit
  real   ,intent(IN):: x1,y1,x2,y2,x3,y3,x4,y4
  write(ggunit,105)x1,y1,x2,y2
  write(ggunit,105)x3,y3,x4,y4
105 format(2(f7.1,' hpu '),'moveto ',2(f7.1,' hpu '),'lineto % ggfmt105') 
end subroutine ggfmt105

subroutine ggfmt106(ggunit,x1,y1,x2,y2)
  integer,intent(IN)::ggunit
  real,intent(IN):: x1,y1,x2,y2
  write(ggunit,106)x1,y1,x2,y2
106 format(2(f7.1,' hpu '),'moveto ',2(f7.1,' hpu '),'lineto % ggfmt106') 
end subroutine ggfmt106

subroutine ggfmt107(ggunit,wid)
  integer,intent(IN):: ggunit
  real,intent(IN):: wid
  write(ggunit,107)wid
107 format(f4.1,' setlinewidth % ggfmt107') 
end subroutine ggfmt107

subroutine ggfmt108(ggunit)
  implicit none
  integer,intent(IN):: ggunit
  write(ggunit,108)
  call gdstroke
108 format('stroke % ggfmt108')
end subroutine ggfmt108

subroutine ggfmt109(ggunit,buf,l)
  implicit none
  integer,intent(IN):: ggunit
  character*80,intent(IN):: buf
  integer,intent(IN):: l
  character*87 buft
  real:: rtod,x,y,xa,ya,a
  data rtod/57.29578/
  call gqdvxy(x,y)
  call gqtang(xa,ya)
  call ggfmt100(ggunit,x,y)
  write(ggunit,*)' % this ggfmt100 called from ggfmt109'
  a=atan2(ya,xa)*rtod
  write(ggunit,109)x,y,a
109 format(f7.1,' hpu ',f7.1,' hpu translate'/f9.3,' rotate % ggfmt109')
  buft(1:1)='('
  buft(2:l+1)=buf(1:l)
  buft(l+2:l+8)=') show%'
  write(ggunit,*)buft(1:l+8)
  write(ggunit,110)a,x,y
110 format('0 0 moveto ',f9.3,' neg rotate'/ &
         f7.1,' hpu neg ',f7.1,' hpu neg translate % ggfmt109/110')
end subroutine ggfmt109

subroutine ggfmt112(y1)
  use pfkinds
  use pfconsts
  use pfstatus
  implicit none
  real,intent(IN):: y1
  integer ifont
  character*40 bufont(nfonts)
  data bufont/            &
       '/Helvetica',      &
       '/Palatino-Roman', &
       '/Courier',        &
       '/Helvetica-Bold', &
       '/Palatino-Bold',  &
       '/Courier-Bold',   &
       '/Symbol'/ 
  write(ggunit,111)bufont(ifonts)
111 format(a40)
  write(ggunit,112)y1
112 format('findfont ',f7.4,' hpc scalefont setfont'/'% ggfmt112')
return
entry glfont
  print '('' The available fonts are:'')'
  do ifont=1,nfonts
     write(6,611)ifont,bufont(ifont)
  enddo
  read(*,*)
611 format(1x,i8,4x,a40)
end subroutine ggfmt112

!subroutine ggfmt125(ggunit,i1,j1,i2,j2)
!  integer,intent(IN):: ggunit,i1,j1,i2,j2
!  write(ggunit,125)i1,i2,j1,j2
!125 format('newpath'/                  &
!         '/left ',i5,' hpu def'/       &
!         '/right ',i5,' hpu def'/      &
!         '/bottom ',i5,' hpu def'/     & 
!         '/top ',i5,' hpu def'/        &
!         'left bottom moveto'/         &
!         'right bottom moveto'/        &
!         'right top moveto'/           &
!         'left top moveto'/            &
!         'closepath clip newpath % ggfmt125') 
!end subroutine ggfmt125

subroutine ggfmt127(ggunit,hue,sat,bri)
  integer,intent(IN):: ggunit
  real,intent(IN):: hue,sat,bri
  write(ggunit,127)hue,sat,bri
127 format(3(1x,f7.5),' sethsbcolor % ggfmt127')
end subroutine ggfmt127

subroutine ggfmt199(ggunit)
  integer,intent(IN):: ggunit
  write(ggunit,199)
199 format('showpage'/'%%EOF % ggfmt199')
end subroutine ggfmt199

subroutine ggjopen(junit)
  common/ggjplot/jrc,j,ibuf(20)
  include 'pathmapb.inc'
  print '('' unit'',i3,'' opened'')',junit
  jrc=1
  j=0
  return
820 print '('' error in ggjopen'')'
end subroutine ggjopen

subroutine gcoast(junit,n,xb,yb,zb,kpen)
  common/ggjplot/jrc,j,ibuf(20)
  dimension xb(*),yb(*),zb(*)
  if(kpen.eq.0)call ggjopen(3)
  kpen=1
  n=0
400 if(j.eq.0)then
     read(junit,rec=jrc,iostat=jstat,err=820) (ibuf(k),k=1,20)
     jrc=jrc+1
  endif
  do i=j+1,20
     n=n+1
     it=ibuf(i)
     it2=it/2
     ipen=it-it2*2
     it4=it2/2
     north=it2-it4*2
     iy=it4/2**14
     ix=it4-iy*2**14
     x=(ix+.5)/2**13-1.
     y=(iy+.5)/2**13-1.
     s=2./(1.+x**2+y**2)
     if(north.eq.1)then
	xb(n)=s*x
	yb(n)=s*y
	zb(n)=s-1.
     else
	xb(n)=-s*x
	yb(n)=y*s
	zb(n)=1.-s
     endif
     if(ipen.eq.0)then
	if(n.eq.1)then
           close(junit)
           print '('' unit '',i3,'' closed'')',junit
           kpen=0
	endif
	j=i
	return
     endif
   enddo
   j=0
   goto 400
820 print '('' error in gcoast'')'; stop
end subroutine gcoast

subroutine ggjwrite(junit,n,dlo,dla)
  common/ggjplot/jrc,j,ibuf(20)
  dimension dlo(*),dla(*)
  pi=4.*atan(1.)
  dtor=pi/180.
  ipen=1
  do i=1,n
     if(i.eq.n)ipen=0
     rlo=dlo(i)*dtor
     rla=dla(i)*dtor
     clo=cos(rlo)
     slo=sin(rlo)
     cla=cos(rla)
     sla=sin(rla)
     x=cla*clo
     y=cla*slo
     z=sla
     if(z.gt.0.)then
	north=1
	si=1./(1.+z)
	x=x*si
	y=y*si
     else
	north=0
	si=1./(1.-z)
	x=-x*si
	y=y*si
     endif
     ix=(x+1.)*2**13
     iy=(y+1.)*2**13
     it=ix+iy*2**14
     it=it*2+north
     it=it*2+ipen
     j=j+1
     if(j>20)then
        write(junit,rec=jrc,iostat=jstat,err=820)(ibuf(k),k=1,20)
        j=1
        jrc=jrc+1
     endif
     ibuf(j)=it
  enddo
  return
820 print '('' error in ggjwrite'')'
end subroutine ggjwrite

subroutine ggjclose(junit)
  common/ggjplot/jrc,j,ibuf(20)
  dimension dlo(1),dla(1)
  dlo(1)=0.
  dla(1)=90.
  call ggjwrite(junit,1,dlo,dla)  ! write a degenerate 1-point line
  do k=j+1,20
     ibuf(k)=0
  enddo
  write(junit,rec=jrc,iostat=jstat,err=820)(ibuf(k),k=1,20)
  close(junit)
  return
820 print '('' error in ggjclose'')'
end subroutine ggjclose
end module pfig