! 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)x1 .or. y(1)y1)z(1)=-1 if(x(-1)x1 .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