program main !*****************************************************************************80 ! !! MAIN is the main program for CRYSTAL_PLOT. ! ! Discussion: ! ! CRYSTAL_PLOT creats plots from the output of the CRYSTAL program. ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! implicit none ! integer, parameter :: maxbot = 20 integer, parameter :: maxl = 64 integer, parameter :: maxobj = 46 ! integer, parameter :: maxm = 2*maxl-1 ! real b1jbl(maxm) real b2jbl(maxm) real b3jbl(maxl) character ( len = 40 ) command real cost real cost2 character ( len = 8 ) date real delx real dely character ( len = 10 ) dev real dxcdp(maxl,maxm) real dxdp(maxl,maxm) real dycdp(maxl,maxm) real dydp(maxl,maxm) real e(maxl,maxm) logical echo character ( len = 80 ) fildat character ( len = 80 ) filgrf character ( len = 80 ) filinp real gamt(maxl,maxm) real grace integer i integer icmax integer icmin integer icolor(maxobj) integer icrys1 integer icrys2 integer ierror integer iplot integer itable integer itemp integer jcmax integer jcmin integer jcrys1 integer jcrys2 integer jtemp integer kcrys(maxl,maxm) integer kmelt(maxl,maxm) integer kvoid(maxl,maxm) integer l logical lbar integer lenc logical lflag(maxl,maxm) logical lgopen logical lnei integer m integer nbot logical ncflag(maxl,maxm) integer ncon logical nflag(maxl,maxm) logical npflag(maxl,maxm) integer nxskip integer nyskip character ( len = 25 ) object(maxobj) logical ovrlay real p(maxl,maxm) real pc(maxl,maxm) real psi(maxl,maxm) logical reflec real rueta(maxl,maxm) real ruksi(maxl,maxm) logical s_eqi real s1(maxl,maxm) real s2(maxl,maxm) real scalecv real scalenb real scalenc real scalenp real scalev logical show(maxobj) real srange real t(maxl,maxm) real te(maxl,maxm) real temp real tempa(maxl,maxm) character ( len = 10 ) time character ( len = 80 ) title character ( len = 80 ) title2 real tk(maxl,maxm) real tnow real tnow2 real u(maxl,maxm) real v(maxl,maxm) logical vpflag(maxl,maxm) real vmag(maxl,maxm) real vort(maxl,maxm) real w(maxl,maxm) real x1max real x1min real x2max real x2min real xbot(maxbot) real xc(maxl,maxm) real xmax real xmin real xp(maxl,maxm) real xsmax real xsmin real xtmax real xtmin real y1max real y1min real y2max real y2min real ybot(maxbot) real yc(maxl,maxm) real ymax real ymin real yp(maxl,maxm) real ysmax real ysmin real ytmax real ytmin ! call timestamp ( ) call hello ! ! Initialize data. ! call init(b1jbl,b2jbl,b3jbl,cost,cost2,delx,dely,dev,echo,fildat,filgrf, & filinp,grace,icmax,icmin,icolor,icrys1,icrys2,iplot,itable,jcmax, & jcmin,jcrys1,jcrys2,l,lbar,m,maxbot,maxl,maxm,maxobj,nbot,ncflag, & ncon,nflag,npflag,nxskip,nyskip,object,ovrlay,p,pc,psi,reflec,scalecv, & scalenb,scalenc,scalenp,scalev,show,title,title2,u,v,vpflag,vmag,vort, & x1max,x1min,x2max,x2min,xbot,xc,xmax,xmin,xp,xsmax,xsmin,xtmax,xtmin, & y1max,y1min,y2max,y2min,ybot,yc,ymax,ymin,yp,ysmax,ysmin,ytmax,ytmin) open ( unit = 17, file = filinp, status = 'replace' ) ! ! Get the next command. ! 10 continue write ( *, * ) ' ' write ( *, * ) '? ("H" for help)' read(*,'(a)',end = 50,err=50)command call flushl(command) write(17,'(a)') trim ( command ) if ( echo ) then write(*,'(a)') trim ( command ) end if if ( command == ' ')go to 10 11 continue ! ! AXis: show axis of symmetry. ! if ( s_eqi ( command(1:2),'ax') ) then show(28) = .not.show(28) if ( show(28) ) then write ( *, '(a)' ) 'The axis of symmetry will be shown.' else write ( *, '(a)' ) 'The axis of symmetry will NOT be shown.' end if ! ! BOundary: show boundary. ! else if ( s_eqi ( command(1:2),'bo') ) then show(1) = .not.show(1) if ( show(1) ) then write ( *, '(a)' ) 'The boundary will be shown.' else write ( *, '(a)' ) 'The boundary will NOT be shown.' end if ! ! BAr: Switch display of color bar. ! else if ( s_eqi ( command(1:2),'ba') ) then lbar = .not.lbar if ( lbar ) then write ( *, '(a)' ) 'The color bar will be shown.' else write ( *, '(a)' ) 'The color bar will NOT be shown.' end if ! ! BH: bottom half. ! else if ( s_eqi ( command,'bh') ) then ysmax = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! BL: bottom left quarter. ! else if ( s_eqi ( command,'bl') ) then xsmax = xsmin+0.5*(xsmax-xsmin) ysmax = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! BC: bottom center quarter. ! else if ( s_eqi ( command,'bc') ) then temp = 0.25*(xsmax-xsmin) xsmin = xsmin+temp xsmax = xsmax-temp ysmax = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! BR: bottom right quarter. ! else if ( s_eqi ( command,'br') ) then xsmin = xsmin+0.5*(xsmax-xsmin) ysmax = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! C: choose colors. ! else if ( s_eqi ( command,'c') ) then write ( *, * ) ' ' write ( *, * ) ' Number Color Name' write ( *, * ) ' ' do i = 1,maxobj write(*,'(1x,i2,2x,i3,2x,a)')i,icolor(i),object(i) end do write ( *, * ) ' ' write ( *, * ) 'Enter an object number, and a color number.' read(*,*,end = 50,err=50)itemp,jtemp write(17,*)itemp,jtemp if ( echo ) then write ( *, * ) itemp,jtemp end if if ( 1 <= itemp .and. itemp <= maxobj ) then icolor(itemp) = jtemp else write ( *, * ) 'Your object number was out of bounds.' end if ! ! CC: choose color contour labels ! ! For some strange reason, in order to make the color table ! active, we have to call NEWFRM! ! else if ( s_eqi ( command(1:2),'cc') ) then if ( dev == ' ' ) then write ( *, * ) ' ' write ( *, * ) 'CPLOT - Warning!' write ( *, * ) ' Please use the "DEV" command to choose' write ( *, * ) ' a device, and THEN the "CC" command.' go to 10 end if if ( s_eqi ( command(1:3),'cc = ') ) then read(command(4:),*,err = 50,end=50)itable else write ( *, * ) ' ' write ( *, * ) 'Built in color tables include:' write ( *, * ) ' ' write ( *, * ) '1 low black to high white.' write ( *, * ) '2 low blue to high yellow.' write ( *, * ) '3 low red, high blue, with bands between.' write ( *, * ) '4 low red, yellow, green, blue, high white.' write ( *, * ) '5 low white, blue, green, yellow, high red' write ( *, * ) '6 low blue to high red.' write ( *, * ) '7 linear scale between 2 user colors.' write ( *, * ) '8 linear scale between N user colors.' write ( *, * ) '9 low white to high black.' write ( *, * ) ' ' write ( *, * ) 'Enter a color table index between 1 and 9,' write ( *, * ) 'or 0 to enter a color table from a file.' read(*,*,end = 50,err=50)itable write(17,*)itable if ( echo ) then write ( *, * ) itable end if end if call gettab(dev,echo,filgrf,grace,icmax,icmin,ierror, & iplot,itable,lgopen,ovrlay) if ( itable == 1 .or. itable == 9 ) then jcmax = 200 jcmin = 32 else jcmax = 255 jcmin = 2 end if write ( *, * ) ' ' write ( *, * ) 'Lowest color used will be JCMIN = ',jcmin write ( *, * ) 'Highest color used will be JCMAX = ',jcmax ! ! CH: center half. ! else if ( s_eqi ( command,'ch') ) then temp = 0.25*(xsmax-xsmin) xsmin = xsmin+temp xsmax = xsmax-temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! CP: show corrected pressure contours. ! else if ( s_eqi ( command,'cp') ) then show(5) = .not.show(5) if ( show(5) ) then write ( *, * ) 'Corrected pressures will be shown.' else write ( *, * ) 'Corrected pressures will NOT be shown.' end if ! ! CPC: show corrected pressure color plots. ! else if ( s_eqi ( command,'cpc') ) then show(12) = .not.show(12) if ( show(12) ) then write ( *, * ) 'Corrected pressure colors will be shown.' else write ( *, * ) 'Corrected pressure colors will NOT be shown.' end if ! ! CRUcible: show crucible wall. ! else if ( s_eqi ( command(1:3),'cru') ) then show(19) = .not.show(19) if ( show(19) ) then write ( *, * ) 'The crucible wall will be shown.' else write ( *, * ) 'The crucible wall will NOT be shown.' end if ! ! CRYB: show the crystal boundary. ! else if ( s_eqi ( command(1:4),'cryb') ) then show(20) = .not.show(20) if ( show(20) ) then write ( *, * ) 'The crystal boundary will be shown.' else write ( *, * ) 'The crystal boundary will NOT be shown.' end if ! ! CRYS: display items in the crystal. ! else if ( s_eqi ( command,'crys') .or. s_eqi ( command,'crystal') ) then show(39) = .not.show(39) if ( show(39) ) then write ( *, * ) 'Items within the crystal will be shown.' else write ( *, * ) 'Items within the crystal will NOT be shown.' end if ! ! CRYSC: display shaded crystal. ! else if ( s_eqi ( command,'crysc') ) then show(44) = .not.show(44) if ( show(44) ) then write ( *, * ) 'The crystal interior will be colored.' else write ( *, * ) 'The crystal interior will NOT be colored.' end if ! ! CV: show control volumes. ! else if ( s_eqi ( command,'cv') ) then show(2) = .not.show(2) if ( show(2) ) then write ( *, * ) 'Control volumes will be shown.' else write ( *, * ) 'Control volumes will NOT be shown.' end if ! ! CVN: show control volume numbers. ! else if ( s_eqi ( command,'cvn') ) then show(43) = .not.show(43) if ( show(43) ) then write ( *, * ) 'Control volume numbers will be shown.' else write ( *, * ) 'Control volume numbers will NOT be shown.' end if ! ! DAT = ! else if ( s_eqi ( command(1:3),'dat') ) then 20 continue if ( s_eqi ( command(1:4),'dat = ') ) then fildat = command(5:) else write ( *, * ) ' ' write ( *, * ) 'Enter the name of the new input data file:' read(*,'(a)',err = 50,end=50)fildat write(17,'(a)') trim ( fildat ) if ( echo ) then write(*,'(a)') trim ( fildat ) end if end if open(unit = 10,file=fildat,status='old',err=60) call rsread(b1jbl,b2jbl,b3jbl,cost,e,gamt,icrys1,icrys2,jcrys1,jcrys2, & kcrys,kmelt,kvoid,l,m,maxbot,maxl,maxm,nbot,p,pc,psi,rueta,ruksi,t, & te,tk,tnow,u,v,vmag,vort,w,xbot,xc,xp,ybot,yc,yp) close(unit = 10) reflec = .false. write ( *, * ) ' ' write ( *, * ) 'CPLOT - Note:' write ( *, * ) ' ' write ( *, * ) ' Region is ',l,' rows wide and ',m,' rows high.' write ( *, * ) ' Time is ',tnow write ( *, * ) ' Cost function is ',cost write ( *, * ) ' ICRYS1 = ',icrys1 write ( *, * ) ' ICRYS2 = ',icrys2 write ( *, * ) ' JCRYS1 = ',jcrys1 write ( *, * ) ' JCRYS2 = ',jcrys2 call rsize(delx,dely,grace,l,m,maxl,maxm,ncflag,npflag,srange,vpflag, & xc,x1max,x1min,x2max,x2min,xmax,xmin,xsmax,xsmin,yc,y1max, & y1min,y2max,y2min,ymax,ymin,ysmax,ysmin) xtmax = xsmax xtmin = xsmin ytmax = ysmax ytmin = ysmin ! ! 'DEV = ' Choose the graphics device. ! else if ( s_eqi ( command(1:3),'dev') ) then if ( dev /= ' ' ) then write ( *, * ) ' ' write ( *, * ) 'CPLOT - Error!' write ( *, * ) ' You have already chosen device '//dev write ( *, * ) ' You may not change your mind!' go to 10 end if if ( s_eqi ( command(1:4),'dev = ') ) then dev = trim ( command(5:) ) else write ( *, * ) ' ' write ( *, * ) 'Enter the graphics device desired.' write ( *, * ) ' ' write ( *, * ) 'Options include:' write ( *, * ) ' ' write ( *, * ) 'CGMB output to a CGM binary file.' write ( *, * ) 'PS output to a PostScript file.' write ( *, * ) 'XWS output to an X window screen.' read(*,'(a)',end = 50,err=50)dev write(17,'(a)') trim ( dev ) if ( echo ) then write(*,'(a)') trim ( dev ) end if end if if ( s_eqi ( dev(1:3),'cgm') ) then dev = 'cgmb' write ( *, * ) 'Output will be to a CGM binary file "cplot.cgm".' else if ( s_eqi ( dev,'ps') ) then write ( *, * ) 'Output will be to a PostScript file "cplot.ps".' else if ( s_eqi ( dev,'xws') ) then write ( *, * ) 'Output will be to an X window screen.' else write ( *, * ) 'Your device '//dev//' was not recognized!' dev = ' ' end if ! ! DIF ! else if ( s_eqi ( command(1:3),'dif') ) then 30 continue write ( *, * ) ' ' write ( *, * ) 'Enter the name of the first input data file:' read(*,'(a)',err = 50,end=50)fildat write(17,'(a)') trim ( fildat ) if ( echo ) then write(*,'(a)') trim ( fildat ) end if open(unit = 10,file=fildat,status='old',err=70) call rsread(b1jbl,b2jbl,b3jbl,cost,e,gamt,icrys1,icrys2,jcrys1,jcrys2, & kcrys,kmelt,kvoid,l,m,maxbot,maxl,maxm,nbot,p,pc,psi,rueta,ruksi,t, & te,tk,tnow,u,v,vmag,vort,w,xbot,xc,xp,ybot,yc,yp) close(unit = 10) reflec = .false. 40 continue write ( *, * ) ' ' write ( *, * ) 'Enter the name of the second input data file:' read(*,'(a)',err = 50,end=50)fildat write(17,'(a)') trim ( fildat ) if ( echo ) then write(*,'(a)') trim ( fildat ) end if open(unit = 10,file=fildat,status='old',err=80) call rsdiff(cost2,dxcdp,dxdp,dycdp,dydp,e,gamt,l,m,maxl,maxm,nbot,p,pc, & psi,rueta,ruksi,t,te,tempa,tk,tnow2,u,v,vmag,vort,w,xc,xp,yc,yp) close(unit = 10) write ( *, * ) 'Region is ',l,' rows wide and ',m,' rows high.' write ( *, * ) 'Time 1 is ',tnow write ( *, * ) 'Time 2 is ',tnow2 write ( *, * ) 'Cost function 1 is ',cost write ( *, * ) 'Cost function 2 is ',cost2 write ( *, * ) 'ICRYS1 = ',icrys1 write ( *, * ) 'JCRYS2 = ',jcrys2 call rsize(delx,dely,grace,l,m,maxl,maxm,ncflag,npflag,srange,vpflag, & xc,x1max,x1min,x2max,x2min,xmax,xmin,xsmax,xsmin,yc,y1max, & y1min,y2max,y2min,ymax,ymin,ysmax,ysmin) ! ! DOUBLE ! else if ( s_eqi ( command(1:3),'dou') ) then reflec = .true. write ( *, * ) ' ' write ( *, * ) 'The region will be reflected around the Y axis.' call double(e,gamt,jcrys1,jcrys2,kcrys,kmelt,kvoid,l,m, & maxbot,maxl,maxm,nbot,ncflag,npflag,p,pc,psi,rueta,ruksi, & t,te,tk,u,v,vmag,vort,vpflag,w,xbot,xc,xp,ybot,yc,yp) write ( *, * ) ' ' write ( *, * ) 'CPLOT - Note:' write ( *, * ) ' The doubled region is ',l,' rows wide and ', m,' rows high.' ! ! Force recomputation of the plot window by setting XSMIN = XSMAX=0. ! xsmin = 0.0 xsmax = 0.0 call rsize(delx,dely,grace,l,m,maxl,maxm,ncflag,npflag, & srange,vpflag,xc,x1max,x1min,x2max,x2min,xmax,xmin,xsmax,xsmin,yc,y1max, & y1min,y2max,y2min,ymax,ymin,ysmax,ysmin) ! ! DXCDP: show (dXCdP,dYCdP) vector plots. ! else if ( s_eqi ( command(1:3),'dxc') ) then show(27) = .not.show(27) if ( show(27) ) then write ( *, * ) '(dXCdP,dYCdP) vector plots will be shown.' else write ( *, * ) '(dXCdP,dYCdP) vector plots will NOT be shown.' end if ! ! ECHO User input will be echoed to the output file. ! else if ( s_eqi ( command,'echo') ) then echo = .not.echo if ( echo ) then write(17,'(a)')command if ( echo ) then write(*,'(a)')command end if write ( *, * ) 'User input will be echoed to the output file.' else write ( *, * ) 'User input will NOT be echoed to the output file.' end if ! ! FILE = : set the name of the graphics output file. ! else if ( s_eqi ( command(1:4),'file') ) then if ( iplot > 0 ) then write ( *, * ) ' ' write ( *, * ) 'Display - Warning:' write ( *, * ) ' It is too late to specify a plot file name.' else if ( s_eqi ( command(1:5),'file = ') ) then filgrf = trim ( command(6:) ) else write ( *, * ) ' ' write ( *, * ) 'Enter the plot file name.' read(*,'(a)',err = 50,end=50)filgrf write(17,'(a)') trim ( filgrf ) if ( echo ) then write(*,'(a)') trim ( filgrf ) end if end if ! ! FRAME: show the frame ! else if ( s_eqi ( command,'frame') ) then show(3) = .true. write ( *, * ) 'A frame will be shown around the picture.' ! ! FS: show free surface. ! else if ( s_eqi ( command(1:2),'fs') ) then show(21) = .not.show(21) if ( show(21) ) then write ( *, * ) 'The free surface will be shown.' else write ( *, * ) 'The free surface will NOT be shown.' end if ! ! FULL: show full picture. ! else if ( s_eqi ( command,'full') ) then ! ! Force recomputation of the plot window by setting XSMIN = XSMAX=0. ! xsmin = 0.0 xsmax = 0.0 call rsize(delx,dely,grace,l,m,maxl,maxm,ncflag,npflag,srange,vpflag, & xc,x1max,x1min,x2max,x2min,xmax,xmin,xsmax,xsmin,yc,y1max, & y1min,y2max,y2min,ymax,ymin,ysmax,ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! GRACE = : set the grace margin. ! else if ( s_eqi ( command(1:5),'grace') ) then if ( s_eqi ( command(1:6),'grace = ') ) then read(command(7:),*,err = 90,end=90)grace else write ( *, * ) 'Enter the grace margin:' read(*,*)grace write(17,*)grace if ( echo ) then write ( *, * ) grace end if end if write ( *, * ) ' ' write ( *, * ) 'CPLOT - Note:' write ( *, * ) ' The grace margin was set to GRACE = ',grace call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! GRAPH ! else if ( s_eqi ( command(1:1),'g') ) then if ( l < 1 .or. m.lt.1 ) then write ( *, * ) ' ' write ( *, * ) 'Whoa! Please read data with the DAT command!' else if ( dev == ' ' ) then write ( *, * ) ' ' write ( *, * ) 'Please use the DEV command first!' else call graph(delx,dely,dev,dxcdp,dycdp,echo,filgrf,gamt,icmax, & icmin,icolor,icrys1,icrys2,iplot,itable,jcmax,jcmin,jcrys1, & jcrys2,kcrys,kmelt,kvoid,l,lbar,lflag,lgopen,m,maxbot,maxl, & maxm,maxobj,nbot,ncflag, & ncon,nflag,npflag,nxskip,nyskip,object,ovrlay,p,pc,psi,s1,s2, & scalecv,scalenb,scalenc,scalenp,scalev,show, & srange,t,te,title,title2,tk,u,v,vpflag,vmag,vort,w, & x1max,x1min,x2max,x2min,xbot,xc,xp,xsmax,xsmin,y1max,y1min, & y2max,y2min,ybot,yc,yp,ysmax,ysmin) end if ! ! HALF ! else if ( s_eqi ( command(1:4),'half') ) then reflec = .false. write ( *, * ) ' ' write ( *, * ) 'The region will NOT be reflected around the Y axis.' call half(e,gamt,jcrys1,jcrys2,kcrys,kmelt,kvoid,l,m, & maxbot,maxl,maxm,nbot,ncflag,npflag,p,pc,psi,rueta,ruksi, & t,te,tk,u,v,vmag,vort,vpflag,w,xbot,xc,xp,ybot,yc,yp) write ( *, * ) ' ' write ( *, * ) 'CPLOT - Note:' write ( *, * ) ' The undoubled region is ',l,' rows wide and ', & m,' rows high.' ! ! Force recomputation of the plot window by setting XSMIN = XSMAX=0. ! xsmin = 0.0 xsmax = 0.0 call rsize(delx,dely,grace,l,m,maxl,maxm,ncflag,nflag,srange,vpflag, & xc,x1max,x1min,x2max,x2min,xmax,xmin,xsmax,xsmin,yc,y1max, & y1min,y2max,y2min,ymax,ymin,ysmax,ysmin) ! ! HELP ! else if ( s_eqi ( command(1:1),'h') ) then call help ! ! ICMAX = : set the maximum available color index. ! else if ( s_eqi ( command(1:6),'icmax = ') ) then read(command(7:),*,err = 50,end=50)icmax if ( icmax > 255 ) then write ( *, * ) 'ICMAX must be no more than 255' icmax = 255 end if write ( *, * ) 'Maximum available color set to ',icmax ! ! ICMIN = : set the minimum available color index. ! else if ( s_eqi ( command(1:6),'icmin = ') ) then read(command(7:),*,err = 50,end=50)icmin if ( icmin < 2 ) then write ( *, * ) 'ICMIN must be no less than 2' icmin = 2 end if write ( *, * ) 'Minimum color set to ',icmin ! ! JCMAX = : set the maximum used color index. ! else if ( s_eqi ( command(1:6),'jcmax = ') ) then read(command(7:),*,err = 50,end=50)jcmax if ( jcmax > 255 ) then write ( *, * ) 'JCMAX must be no more than 255.' jcmax = 255 end if write ( *, * ) 'Maximum used color set to ',jcmax ! ! JCMIN = : set the minimum used color index. ! else if ( s_eqi ( command(1:6),'jcmin = ') ) then read(command(7:),*,err = 50,end=50)jcmin if ( jcmin < 2 ) then write ( *, * ) 'JCMIN must be no less than 2.' jcmin = 2 end if write ( *, * ) 'Minimum used color set to ',jcmin ! ! LH: left half. ! else if ( s_eqi ( command,'lh') ) then xsmax = xsmin+0.5*(xsmax-xsmin) temp = 0.25*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! MC: middle center quarter. ! else if ( s_eqi ( command,'mc') ) then temp = 0.25*(xsmax-xsmin) xsmin = xsmin+temp xsmax = xsmax-temp temp = 0.25*(ysmax-ysmin) ysmax = ysmax-temp ysmin = ysmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! MELT: display items in the melt. ! else if ( s_eqi ( command,'melt') ) then show(40) = .not.show(40) if ( show(40) ) then write ( *, * ) 'Items within the melt will be shown.' else write ( *, * ) 'Items within the melt will NOT be shown.' end if ! ! MELTC: display shaded melt. ! else if ( s_eqi ( command,'meltc') ) then show(45) = .not.show(45) if ( show(45) ) then write ( *, * ) 'The interior of the melt will be colored.' else write ( *, * ) 'The interior of the melt will NOT be colored.' end if ! ! MH: middle half. ! else if ( s_eqi ( command,'mh') ) then temp = 0.25*(ysmax-ysmin) ysmax = ysmax-temp ysmin = ysmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! ML: middle left quarter. ! else if ( s_eqi ( command,'ml') ) then xsmax = xsmin+0.5*(xsmax-xsmin) temp = 0.25*(ysmax-ysmin) ysmax = ysmax-temp ysmin = ysmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! MR: middle right quarter. ! else if ( s_eqi ( command,'mr') ) then xsmin = xsmin+0.5*(xsmax-xsmin) temp = 0.25*(ysmax-ysmin) ysmax = ysmax-temp ysmin = ysmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! NB: show nodes on bottom that define crucible shape. ! else if ( s_eqi ( command,'nb') ) then show(34) = .not.show(34) if ( show(34) ) then write ( *, * ) 'Nodes on bottom will be shown.' else write ( *, * ) 'Nodes on bottom will NOT be shown.' end if ! ! NC: show corner nodes ! else if ( s_eqi ( command,'nc') ) then show(29) = .not.show(29) if ( show(29) ) then write ( *, * ) 'Corner nodes will be shown.' else write ( *, * ) 'Corner nodes will NOT be shown.' end if ! ! NCON = : Set the number of contour lines. ! else if ( s_eqi ( command(1:4),'ncon') ) then if ( s_eqi ( command(1:5),'ncon = ') ) then read(command(6:),*,err = 50,end=50)ncon write ( *, * ) 'Number of contour lines set to ',ncon else write ( *, * ) 'Enter number of contour lines.' read(*,*)ncon write(17,*)ncon if ( echo ) then write ( *, * ) ncon end if end if ! ! NOFRAME: no frame, please ! else if ( s_eqi ( command,'noframe') ) then show(3) = .false. write ( *, * ) 'No frame will be shown around the picture.' ! ! NP: show primary nodes ! else if ( s_eqi ( command,'np') ) then show(4) = .not.show(4) if ( show(4) ) then write ( *, * ) 'Primary nodes will be shown.' else write ( *, * ) 'Primary nodes will NOT be shown.' end if ! ! NXSKIP = : skip value for column nodes. ! else if ( s_eqi ( command(1:7),'nxskip = ') ) then read(command(8:),*,err = 50,end=50)itemp if ( itemp <= 0 ) then write ( *, * ) 'Error! NXSKIP must be 1 or more!' else nxskip = itemp write ( *, * ) 'Skip value for column nodes is ',nxskip end if ! ! NYSKIP = : skip value for row nodes. ! else if ( s_eqi ( command(1:7),'nyskip = ') ) then read(command(8:),*,err = 50,end=50)itemp if ( itemp <= 0 ) then write ( *, * ) 'Error! NYSKIP must be 1 or more!' else nyskip = itemp write ( *, * ) 'Skip value for row nodes is ',nyskip end if ! ! OVERLAY: Switch the overlay value. ! else if ( s_eqi ( command,'overlay') ) then ovrlay = .not.ovrlay if ( ovrlay ) then write ( *, * ) 'Plots will be overlayed until next OVERLAY.' else write ( *, * ) 'This overlay plot is done.' call newfrm end if ! ! P: show pressure contours. ! else if ( s_eqi ( command,'p') ) then show(5) = .not.show(5) if ( show(5) ) then write ( *, * ) 'Pressures will be shown.' else write ( *, * ) 'Pressures will NOT be shown.' end if ! ! PC: show pressure color plots. ! else if ( s_eqi ( command,'pc') ) then show(12) = .not.show(12) if ( show(12) ) then write ( *, * ) 'Pressure colors will be shown.' else write ( *, * ) 'Pressure colors will NOT be shown.' end if ! ! Q = Quit. ! QY = QUIT NOW! ! else if ( s_eqi ( command(1:1),'q') ) then if ( s_eqi ( command(2:2),'y') ) then command = 'y' else write ( *, * ) 'Enter "y" to confirm you want to quit!' read(*,'(a)')command call flushl(command) write(17,'(a)') trim ( command ) if ( echo ) then write(*,'(a)') trim ( command ) end if end if if ( s_eqi ( command,'y') ) then if ( lgopen ) then call grfcls if ( lnei ( dev, 'cgmb' ) ) then call delete ( 'cgmout' ) end if end if close(unit = 17) stop end if ! ! RH: right half. ! else if ( s_eqi ( command,'rh') ) then temp = 0.50*(xsmax-xsmin) xsmin = xsmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! S: show stream lines ! else if ( s_eqi ( command,'s') ) then show(6) = .not.show(6) if ( show(6) ) then write ( *, * ) 'Stream lines will be shown.' else write ( *, * ) 'Stream lines will NOT be shown.' end if ! ! SC: show stream colors ! else if ( s_eqi ( command,'sc') ) then show(35) = .not.show(35) if ( show(35) ) then write ( *, * ) 'Stream colors will be shown.' else write ( *, * ) 'Stream colors will NOT be shown.' end if ! ! SCALENB = : set the bottom node scale. ! else if ( s_eqi ( command(1:8),'scalenb = ') ) then read(command(9:),*,err = 50,end=50)scalenb write ( *, * ) 'Scale factor for bottom nodes set to ',scalenb ! ! SCALENC = : set the corner node scale. ! else if ( s_eqi ( command(1:8),'scalenc = ') ) then read(command(9:),*,err = 50,end=50)scalenc write ( *, * ) 'Scale factor for corner nodes set to ',scalenc ! ! SCALENP = : set the primary node scale. ! else if ( s_eqi ( command(1:8),'scalenp = ') ) then read(command(9:),*,err = 50,end=50)scalenp write ( *, * ) 'Scale factor for primary nodes set to ',scalenp ! ! SCALEV = : set the velocity vector scale. ! else if ( s_eqi ( command(1:7),'scalev = ') ) then read(command(8:),*,err = 50,end=50)scalev write ( *, * ) 'Velocity vector scale factor set to ',scalev ! ! TC: top center quarter. ! else if ( s_eqi ( command,'tc') ) then temp = 0.25*(xsmax-xsmin) xsmin = xsmin+temp xsmax = xsmax-temp ysmin = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! TE: show turbulent epsilon contours. ! else if ( s_eqi ( command,'te') ) then show(22) = .not.show(22) if ( show(22) ) then write ( *, * ) 'Turbulent epsilon will be shown.' else write ( *, * ) 'Turbulent epsilon will NOT be shown.' end if ! ! TEMP: show temperature contours. ! else if ( s_eqi ( command,'temp') ) then show(11) = .not.show(11) if ( show(11) ) then write ( *, * ) 'Temperatures will be shown.' else write ( *, * ) 'Temperatures will NOT be shown.' end if ! ! TEMPC: show temperature colors. ! else if ( s_eqi ( command,'tempc') ) then show(13) = .not.show(13) if ( show(13) ) then write ( *, * ) 'Temperature colors will be shown.' else write ( *, * ) 'Temperature colors will NOT be shown.' end if ! ! TH: top half ! else if ( s_eqi ( command,'th') ) then ysmin = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! TITLE2 = : enter subtitle ! else if ( s_eqi ( command(1:7),'title2 = ') ) then title2 = trim ( command(8:) ) else if ( s_eqi ( command(1:6),'title2') ) then write ( *, * ) 'Enter the subtitle:' read(*,'(a)',end = 50)title2 write(17,'(a)') trim ( title2 ) if ( echo ) then write(*,'(a)') trim ( title2 ) end if ! ! TITLE = : enter title ! else if ( s_eqi ( command(1:6),'title = ') ) then title = trim ( command(7:) ) else if ( s_eqi ( command(1:5),'title') ) then write ( *, * ) 'Enter the plot title:' read(*,'(a)',end = 50)title write(17,'(a)') trim ( title ) if ( echo ) then write(*,'(a)') trim ( title ) end if ! ! TK: show turbulent kinetic energy contours. ! else if ( s_eqi ( command,'tk') ) then show(23) = .not.show(23) if ( show(23) ) then write ( *, * ) 'Turbulent KE will be shown.' else write ( *, * ) 'Turbulent KE will NOT be shown.' end if ! ! TL: top left quarter. ! else if ( s_eqi ( command,'tl') ) then xsmax = xsmin+0.5*(xsmax-xsmin) ysmin = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! TR: top right quarter. ! else if ( s_eqi ( command,'tr') ) then xsmin = xsmin+0.5*(xsmax-xsmin) ysmin = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! U: show unit velocity direction vectors. ! else if ( s_eqi ( command,'u') ) then show(9) = .not.show(9) if ( show(9) ) then write ( *, * ) 'Unit velocity direction field will be shown.' else write ( *, * ) 'Unit velocity direction field will NOT be shown.' end if ! ! V: show velocity vectors. ! else if ( s_eqi ( command,'v') ) then show(8) = .not.show(8) if ( show(8) ) then write ( *, * ) 'Velocities will be shown.' else write ( *, * ) 'Velocities will NOT be shown.' end if ! ! VIS: show viscosity contours. ! else if ( s_eqi ( command,'vis') ) then show(24) = .not.show(24) if ( show(24) ) then write ( *, * ) 'Viscosity contours will be in next plot.' else write ( *, * ) 'Viscosity contours will NOT be in next plot.' end if ! ! VM: show velocity magnitude contours. ! else if ( s_eqi ( command,'vm') ) then show(36) = .not.show(36) if ( show(36) ) then write ( *, * ) 'Velocity magnitudes will be shown.' else write ( *, * ) 'Velocity magnitudes will NOT be shown.' end if ! ! VMC: show velocity magnitude colors. ! else if ( s_eqi ( command,'vmc') ) then show(37) = .not.show(37) if ( show(37) ) then write ( *, * ) 'Velocity magnitude colors will be shown.' else write ( *, * ) 'Velocity magnitude colors will NOT be shown.' end if ! ! VOID: display items in the void. ! else if ( s_eqi ( command,'void') ) then show(38) = .not.show(38) if ( show(38) ) then write ( *, * ) 'Items in the void will be shown.' else write ( *, * ) 'Items in the void will NOT be shown.' end if ! ! VOIDC: display shaded void. ! else if ( s_eqi ( command,'voidc') ) then show(46) = .not.show(46) if ( show(46) ) then write ( *, * ) 'The void will be colored.' else write ( *, * ) 'The void will NOT be colored.' end if ! ! VORT: show vorticity contours. ! else if ( s_eqi ( command,'vort') ) then show(41) = .not.show(41) if ( show(41) ) then write ( *, * ) 'Vorticity contours will be shown.' else write ( *, * ) 'Vorticity contours will NOT be shown.' end if ! ! VORTC: show vorticity colors. ! else if ( s_eqi ( command,'vortc') ) then show(42) = .not.show(42) if ( show(42) ) then write ( *, * ) 'Vorticity colors will be shown.' else write ( *, * ) 'Vorticity colors will NOT be shown.' end if ! ! VPN: set visible primary nodes. ! else if ( s_eqi ( command,'vpn') ) then call vizpn(echo,l,m,maxl,maxm,vpflag,xp,yp) ! ! W: show angular momentum contours. ! else if ( s_eqi ( command,'w') ) then show(10) = .not.show(10) if ( show(10) ) then write ( *, * ) 'Angular momentum contours will be in next plot.' else write ( *, * ) 'Angular momentum contours will NOT be shown.' end if ! ! WC: show angular momentum color plots. ! else if ( s_eqi ( command,'wc') ) then show(14) = .not.show(14) if ( show(14) ) then write ( *, * ) 'Angular momentum colors will be shown.' else write ( *, * ) 'Angular momentum colors will NOT be shown.' end if ! ! X: set the data window. ! else if ( s_eqi ( command,'x') ) then call getwin(echo,grace,srange,xmax,xmin,x1max,x1min,x2max, & x2min,xsmax,xsmin,ymax,ymin,y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! XC: show X coordinate contours. ! else if ( s_eqi ( command,'xc') ) then show(30) = .not.show(30) if ( show(30) ) then write ( *, * ) 'X coordinate contours will be shown.' else write ( *, * ) 'X coordinate contours will NOT be shown.' end if ! ! XCC: show X coordinate colors. ! else if ( s_eqi ( command,'xcc') ) then show(31) = .not.show(31) if ( show(31) ) then write ( *, * ) 'X coordinate colors will be shown.' else write ( *, * ) 'X coordinate colors will NOT be shown.' end if ! ! YC: show Y coordinate contours. ! else if ( s_eqi ( command,'yc') ) then show(32) = .not.show(32) if ( show(32) ) then write ( *, * ) 'Y coordinate contours will be shown.' else write ( *, * ) 'Y coordinate contours will NOT be shown.' end if ! ! YCC: show Y coordinate colors. ! else if ( s_eqi ( command,'ycc') ) then show(33) = .not.show(33) if ( show(33) ) then write ( *, * ) 'Y coordinate colors will be shown.' else write ( *, * ) 'Y coordinate colors will NOT be shown.' end if ! ! Unrecognized command. ! else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CPLOT - Warning!' write ( *, '(a)' ) ' Unrecognized command:' // trim ( command ) end if go to 10 50 continue write ( *, * ) ' ' write ( *, * ) 'CPLOT - Fatal error!' write ( *, * ) ' Input error or end of file.' write ( *, * ) ' ' write ( *, * ) ' Warning! The graphics file may not be' write ( *, * ) ' properly terminated.' stop 60 continue write ( *, * ) ' ' write ( *, * ) 'CPLOT - Warning!' write ( *, * ) ' The file you named could not be opened!' fildat = ' ' command = 'dat' go to 11 70 continue write ( *, * ) ' ' write ( *, * ) 'CPLOT - Warning!' write ( *, * ) ' The file you named could not be opened!' fildat = ' ' command = 'dif' go to 11 80 continue write ( *, * ) ' ' write ( *, * ) 'CPLOT - Warning!' write ( *, * ) ' The file you named could not be opened!' fildat = ' ' command = 'dif' go to 11 90 continue write ( *, * ) ' ' write ( *, * ) 'CPLOT - Warning!' write ( *, * ) ' Your input was not understood.' go to 10 end subroutine arrow ( xstart, ystart, xtip, ytip ) ! !*****************************************************************************80 ! !! ARROW can be used to draw an arrow at any point on a graph. ! ! ! Discussion: ! ! The arrow will stretch between two user specified points. ! ! The "head" of the arrow may be fatter or thinner than expected ! if the X and Y scales of the graph are not in the same ! proportions. ! ! left ! |\ ! | \ ! start ------- base tip ! | / ! |/ ! rite ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real XSTART, YSTART, the starting point for the ! arrow. ! ! Input, real XTIP, YTIP, the end point for the arrow. ! implicit none ! real alpha real del real dist real, parameter :: pi = 3.1415926E+00 real theta real xbase real xleft real xrite real xstart real xtip real ybase real yleft real yrite real ystart real ytip ! if ( xstart == xtip .and. ystart == ytip)return theta = 0.5*pi-atan2(2.0,1.0) dist = sqrt((xtip-xstart)**2+(ytip-ystart)**2) alpha = atan2(ytip-ystart,xtip-xstart) del = sqrt(5.0)*dist/3.0 call movcgm(xstart,ystart) xbase = xstart+dist*cos(alpha)*2.0/3.0 ybase = ystart+dist*sin(alpha)*2.0/3.0 call drwcgm(xbase,ybase) xleft = xstart+del*cos(alpha-theta) yleft = ystart+del*sin(alpha-theta) call drwcgm(xleft,yleft) call drwcgm(xtip,ytip) xrite = xstart+del*cos(alpha+theta) yrite = ystart+del*sin(alpha+theta) call drwcgm(xrite,yrite) call drwcgm(xbase,ybase) return end subroutine box(xmin,xmax,ymin,ymax) ! !*****************************************************************************80 ! !! BOX draws a rectangle whose corners are specified by the user. ! ! ! Discussion: ! ! The rectangle drawn by box has the corners: ! ! (XMIN,YMAX) (XMAX,YMAX) ! ! (XMIN,YMIN) (XMAX,YMIN) ! ! BOX can be used to place a rectangle anywhere in the picture. ! However, BOX may also be used to place a rectangle around the ! entire picture, producing a "frame". ! ! The DRAWCGM routine PLYLIN is used to draw the box, and hence ! the color of the line may be set by calling the DRAWCGM routine ! LINCLR first. ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real XMIN, XMAX, the minimum and maximum X ! coordinates of the box. ! ! Input, real YMIN, YMAX, the minimum and maximum Y ! coordinates of the box. ! implicit none ! integer, parameter :: npoints = 5 ! real x(npoints) real xmax real xmin real y(npoints) real ymax real ymin ! x(1) = xmin y(1) = ymin x(2) = xmax y(2) = ymin x(3) = xmax y(3) = ymax x(4) = xmin y(4) = ymax x(5) = xmin y(5) = ymin call plylin(npoints,x,y) return end subroutine buzz(dev,x1max,x1min,y1max,y1min) ! !*****************************************************************************80 ! !! BUZZ forces the graphical system to slow down a bit. ! ! ! Discussion: ! ! BUZZ is just "busy work" which seems to fix a problem that occurs ! when the XWS interface is used. In that case, the last bit of the ! graph is not drawn. So here, we just make the last bit of the graph ! something we don't care about. ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character*10 DEV, the graphics output device to be used. ! Current legal include: ! ! cgmb - CGM binary file. ! ps - PostScript file. ! xws - X window screen (interactive). ! ! Output, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! Output, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! implicit none ! character ( len = 10 ) dev integer i integer icolor real x1max real x1min real y1max real y1min ! if ( dev == 'xws' ) then icolor = 0 call linclr(icolor) do i = 1,100 call box(x1min,x1max,y1min,y1max) end do end if return end subroutine cbar(icolor,jcmax,jcmin,maxobj,ncon,smax,smin,srange,x1,x2,y1,y2) ! !*****************************************************************************80 ! !! CBAR draws a color bar. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! ICOLOR Input, integer ICOLOR(MAXOBJ). ! Contains the color indexes for each object. ! However, in some cases, ICOLOR is actual a color table ! index. ! ! JCMAX, ! JCMIN Input, integer JCMAX, JCMIN, the maximum and ! minimum color indices to use in the color bar. ! ! MAXOBJ Input, integer MAXOBJ. ! The number of graphical "objects". ! ! NCON Input, integer NCON, the number of color contour ! regions drawn, and hence, the number of colors ! to be displayed in the color bar. ! ! SMAX, ! SMIN Input, real SMAX, SMIN, the maximum and minimum ! values of the quantity whose color contours are ! being drawn. These numbers will be printed along ! with the color bar. ! ! SRANGE Input, real SRANGE. ! The maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! X1, ! X2, ! Y1, ! Y2 Input, real X1, X2, Y1, Y2, specify the minimum and ! maximum X and Y coordinates of the color bar. ! implicit none ! integer maxobj ! real angle character chrrel*14 character chrtmp*14 real cwide character ( len = 6 ) flush integer i integer icolor(maxobj) integer jcmax integer jcmin integer jcolor integer lent integer ncon real pwide real smax real smin real srange real x real x1 real x2 real xcorn(5) real xl real xr real y real y1 real y2 real ycorn(5) ! ycorn(1) = y1 ycorn(2) = y1 ycorn(3) = y2 ycorn(4) = y2 ycorn(5) = y1 call linclr(icolor(1)) do i = 0,ncon xl = ((ncon+1-i)*x1+i*x2)/real(ncon+1) xr = ((ncon-i)*x1+(i+1)*x2)/real(ncon+1) xcorn(1) = xl xcorn(2) = xr xcorn(3) = xr xcorn(4) = xl xcorn(5) = xl jcolor = int(((ncon-i)*jcmin+i*jcmax)/real(ncon)) call filclr(jcolor) call plygon(4,xcorn,ycorn) call plylin(5,xcorn,ycorn) end do ! ! Print labels for the lowest and highest contours. ! cwide = 0.9*srange/40.0 chrtmp = chrrel(smin) call s_blank_delete ( chrtmp ) angle = 0.0 x = x1 y = y1-1.5*cwide pwide = srange flush = 'left' call s_plot(angle,cwide,pwide,chrtmp,x,y,flush) chrtmp = chrrel(smax) call s_blank_delete ( chrtmp ) angle = 0.0 x = x2 y = y1-1.5*cwide pwide = srange flush = 'right' call s_plot(angle,cwide,pwide,chrtmp,x,y,flush) return end subroutine cbox(grace,x1max,x1min,y1max,y1min) ! !*****************************************************************************80 ! !! CBOX draws a 16 by 16 color box. ! ! ! Discussion: ! ! You may want to call NEWFRM before or after calling CBOX, in ! order to clear the frame! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! GRACE Input, real GRACE. ! The size of the "grace" margin on the plot. ! ! X1MAX, ! X1MIN Output, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! Y1MAX, ! Y1MIN Output, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! implicit none ! integer, parameter :: npoly = 5 ! real grace integer i integer icolor integer ierror integer j real x(npoly) real x1max real x1min real x2max real x2min real y(npoly) real y1max real y1min real y2max real y2min ! ! Set the coordinate system to be 0 < = X <= 16.0, ! and 0.0 < = Y <= 16.0 ! x2min = 0.0E+00 x2max = 16.0E+00 y2min = 0.0E+00 y2max = 16.0E+00 x1min = x2min-grace*(x2max-x2min) x1max = x2max+grace*(x2max-x2min) y1min = y2min-grace*(y2max-y2min) y1max = y2max+grace*(y2max-y2min) call setwcd(x1min,y1min,x1max,y1max,ierror) ! ! Draw the color boxes. ! icolor = 0 do i = 1,16 y(1) = 16-i y(2) = 16-i y(3) = 17-i y(4) = 17-i do j = 1,16 call filclr(icolor) icolor = icolor+1 x(1) = j-1 x(2) = j x(3) = j x(4) = j-1 call plygon(4,x,y) end do end do ! ! Draw black lines around the boxes. ! icolor = 1 call linclr(icolor) do i = 1,16 y(1) = 16-i y(2) = 17-i y(3) = 17-i y(4) = 16-i y(5) = y(1) do j = 1,16 x(1) = j-1 x(2) = j-1 x(3) = j x(4) = j x(5) = x(1) call plylin(5,x,y) end do end do return end subroutine ch_cap ( c ) ! !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none ! character c integer itemp ! itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end function chrrel(rval) ! !*********************************************************************** ! !! CHRREL accepts a real number in RVAL and returns in CHRREL a ! 14-character right-justified representation of that number. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! RVAL Input, real RVAL, a real number. ! ! CHRREL Output (through function value), character*14 CHRREL, ! a right-justified character variable containing the ! representation of RVAL, using a G14.6 format. ! implicit none ! character ( len = 14 ) chrrel character ( len = 14 ) chrtmp real rval ! ! We can't seem to write directly into CHRREL because of compiler ! quibbles. ! if ( real(int(rval)) == rval .and. abs(rval) < 1.0E+13 ) then write(chrtmp,'(i14)')int(rval) else write(chrtmp,'(g14.6)')rval end if chrrel = chrtmp return end subroutine colcon(a,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) ! !*****************************************************************************80 ! !! COLCON supervises the creation of a color contour plot. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! A Input, real A(MAXL,MAXM), the quantity whose contour ! is desired. ! ! ANAME Input, character*(*) ANAME, the name of the quantity ! to be contoured. ! ! ICOLOR Input, integer ICOLOR(MAXOBJ). ! Contains the color indexes for each object. ! However, in some cases, ICOLOR is actual a color table ! index. ! ! JCMAX, ! JCMIN Input, integer JCMAX, JCMIN, the maximum and ! minimum color indices to use in the color bar. ! ! L Input, integer L, the number of rows of data. ! ! LBAR Input, logical LBAR, is .TRUE. if the color bar should ! be shown. ! ! M Input, integer M, the number of columns of data. ! ! MAXL, ! MAXM Input, integer MAXL, MAXM, the maximum allowed number of ! rows and columns of data. ! ! MAXOBJ Input, integer MAXOBJ. ! The number of graphical "objects". ! ! NCON Input, integer NCON. ! The number of contour lines to be drawn. ! ! NFLAG Input, logical NFLAG(MAXL,MAXM). ! ! NFLAG is used to "flag" which nodes are active, ! that is, to be displayed, and which not, in the graph. ! ! SRANGE Input, real SRANGE. ! The maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! X Input, real X(MAXL,MAXM). ! The X coordinates of the nodes. ! ! X2MAX, ! X2MIN Input, real X2MAX, X2MIN, the maximum and minimum X ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! Y Input, real Y(MAXL,MAXM). ! The Y coordinates of the nodes. ! ! Y1MAX Input, real Y1MAX, the maximum Y coordinates of the plot, ! which includes a small grace margin. ! ! Y2MAX Input, real Y2MAX, the maximum Y coordinates that should be ! used for plotting. ! implicit none ! integer maxl integer maxm integer maxobj ! real a(maxl,maxm) character ( len = * ) aname logical echo integer icolor(maxobj) integer jcmax integer jcmin integer l logical lbar integer m integer ncon logical nflag(maxl,maxm) real smax real smin real srange real xp(maxl,maxm) real x1 real x2 real x2max real x2min real yp(maxl,maxm) real y1 real y1max real y2 real y2max ! ! Get the minimum and maximum values of the data. ! call fsize(l,m,maxl,maxm,nflag,a,smax,smin) ! ! If there is no variation in the data, don't try to draw a ! contour plot. ! if ( smax <= smin ) then write ( *, * ) ' ' write ( *, * ) 'COLCON - Warning!' write ( *, * ) ' No color contours, all values equal.' return end if ! ! Allow the user to adjust the range of the data. ! call setsiz(echo,aname,smax,smin) ! ! Draw the color contour plot. ! call tricol(jcmax,jcmin,l,m,maxl,maxm,ncon,nflag,a,smax,smin,xp,yp) ! ! Draw the color bar. ! if ( lbar ) then x1 = x2min x2 = x2max y1 = y2max y2 = y1max call cbar(icolor,jcmax,jcmin,maxobj,ncon,smax,smin,srange,x1,x2,y1,y2) end if return end subroutine cross(px,py,qx,qy,sl,sm,sh,sval,xl,xm,xh,yl,ym,yh) ! !*****************************************************************************80 ! !! CROSS finds two places where a given value occurs on a triangle. ! ! ! Discussion: ! ! The corners of the triangle are (XL,YL), (XM,YM) and ! (XH,YH), and the associated S values are SL, SM and SH. It ! must be the case that SL < = SM <= SH. ! ! CROSS returns two points: ! ! (PX,PY), which occurs on one of the two sides that include ! (XM,YM), and ! (QX,QY), which occurs on the side between (XL,YL) and (XH,YH). ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! PX, ! PY Output, real PX, PY, the X and Y coordinates of a point ! at which S = SVAL, lying on a side of the triangle which ! ends at (XM,YM). ! ! QX, ! QY Output, real QX, QY, the X and Y coordinates of a point ! at which S = SVAL, lying on the side of the triangle which ! lies between (XL,YL) and (XH,YH). ! ! SL, ! SM, ! SH Input, real SL, SM, SH, the low, medium, and high values ! of S, associated with the three corners. ! ! SVAL Input, real SVAL, the value of S for which a contour line ! is sought. ! ! XL, ! XM, ! XH Input, real XL, XM, XH, the X coordinates of the nodes ! at which the low, medium and high values of S occur. ! ! YL, ! YM, ! YH Input, real YL, YM, YH, the Y coordinates of the nodes ! at which the low, medium and high values of S occur. ! implicit none ! real px real py real qx real qy real sl real sm real sh real sval real xl real xm real xh real yl real ym real yh ! if ( sval < sl ) then px = xl py = yl qx = xl qy = yl else if ( sval >= sh ) then px = xh py = yh qx = xh qy = yh else if ( sval < sm ) then px = xl+(sval-sl)*(xm-xl)/(sm-sl) py = yl+(sval-sl)*(ym-yl)/(sm-sl) else px = xm+(sval-sm)*(xh-xm)/(sh-sm) py = ym+(sval-sm)*(yh-ym)/(sh-sm) end if qx = xl+(sval-sl)*(xh-xl)/(sh-sl) qy = yl+(sval-sl)*(yh-yl)/(sh-sl) end if return end subroutine delete(filnam) ! !*****************************************************************************80 ! !! DELETE deletes a file. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character*(*) FILNAM, the file to be deleted. ! implicit none ! character ( len = * ) filnam ! open(unit = 99,file=filnam,status='old',err=10) close(unit = 99,status='delete',err=10) 10 continue return end subroutine diamnd(xcentr,ycentr,radius,filled) ! !*****************************************************************************80 ! !! DIAMND may be used to draw an open or filled diamond of a given size. ! ! ! Discussion: ! ! DIAMND calls PLYLIN to draw an open diamond, or PLYGON to draw a ! filled diamond. ! ! To control the color of the diamond, simply call the DRAWCGM ! routine LINCLR before drawing open diamonds or FILCOR before ! drawing closed diamonds. ! ! If the X and Y coordinate systems do not have the same scale, ! the diamond drawn will be "flattened". This can happen if the ! routine SETWCD has been used to set the X and Y dimensions to ! different extents. ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real XCENTR, YCENTR, the X and Y coordinates of ! the center of the diamond. ! ! Input, real RADIUS, the radius of the diamond. ! ! Input, logical FILLED, .TRUE. if the diamond is to be filled. ! implicit none ! integer, parameter :: npts = 5 ! logical filled real radius real xcentr real xpoint(npts) real ycentr real ypoint(npts) ! xpoint(1) = xcentr+radius ypoint(1) = ycentr xpoint(2) = xcentr ypoint(2) = ycentr+radius xpoint(3) = xcentr-radius ypoint(3) = ycentr xpoint(4) = xcentr ypoint(4) = ycentr-radius xpoint(5) = xcentr+radius ypoint(5) = ycentr if ( filled ) then call plygon(npts-1,xpoint,ypoint) else call plylin(npts,xpoint,ypoint) end if return end subroutine double(e,gamt,jcrys1,jcrys2,kcrys,kmelt,kvoid,l,m, & maxbot,maxl,maxm,nbot,ncflag,npflag,p,pc,psi,rueta,ruksi, & t,te,tk,u,v,vmag,vort,vpflag,w,xbot,xc,xp,ybot,yc,yp) ! !*****************************************************************************80 ! !! DOUBLE "reflects" the data around the Y axis. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real E(MAXL,MAXM), the magnetic stream function. ! ! Input/output, real GAMT(MAXL,MAXM), the diffusion coefficient. ! ! Input/output, integer JCRYS1, the J coordinate of the column of ! corner nodes that define the left side of the crystal. ! ! Input/output, integer JCRYS2, the J coordinate of the column of ! corner nodes that define the right side of the crystal. ! ! Input/output, integer KCRYS(MAXL,MAXM), ! 0, if control volume (I,J) is away from the crystal. ! 1, if control volume (I,J) is on the external crystal boundary. ! 2, if control volume (I,J) is on the internal crystal boundary. ! 3, if control volume (I,J) is in the crystal interior. ! ! Input/output, integer KMELT(MAXL,MAXM), ! 0, if control volume (I,J) is away from the melt. ! 1, if control volume (I,J) is on the external melt boundary. ! 2, if control volume (I,J) is on the internal melt boundary. ! 3, if control volume (I,J) is in the melt interior. ! ! Input/output, integer KVOID(MAXL,MAXM), ! 0, if control volume (I,J) is away from the void. ! 1, if control volume (I,J) is on the external void boundary. ! 2, if control volume (I,J) is on the internal void boundary. ! 3, if control volume (I,J) is in the void interior. ! ! Input, integer L, the number of rows of data. ! ! Input/output, integer M, the number of columns of data. ! ! Input, integer MAXL, MAXM, the maximum allowed number of ! rows and columns of data. ! ! Input/output, real P(MAXL,MAXM), the pressure. ! ! Input/output, real PC(MAXL,MAXM), the corrected pressure. ! ! Input/output, real PSI(MAXL,MAXM), the stream function. ! ! Input/output, real RUETA(MAXL,MAXM), RUKSI(MAXL,MAXM), the ! the momentum in the ETA and KSI directions. ! ! Input/output, real T(MAXL,MAXM), the temperature. ! ! Input/output, real TE(MAXL,MAXM), the turbulent epsilon. ! ! Input/output, real TK(MAXL,MAXM), the turbulent K. ! ! Input/output, real U(MAXL,MAXM), the horizontal velocity. ! ! Input/output, real V(MAXL,MAXM), the vertical velocity. ! ! Input/output, real W(MAXL,MAXM), the axial velocity. ! ! Input/output, real X(MAXL,MAXM), the X coordinates of the primary ! nodes. ! ! Input/output, real XC(MAXL,MAXM), the X coordinates of the corner ! nodes. ! ! Input/output, real Y(MAXL,MAXM), the Y coordinates of the primary ! nodes. ! ! Input/output, real YC(MAXL,MAXM), the Y coordinates of the corner ! nodes. ! implicit none ! integer maxbot integer maxl integer maxm ! real e(maxl,maxm) real gamt(maxl,maxm) integer i integer j integer jcrys1 integer jcrys2 integer kcrys(maxl,maxm) integer kmelt(maxl,maxm) integer kvoid(maxl,maxm) integer l integer m integer nbot logical ncflag(maxl,maxm) logical npflag(maxl,maxm) real p(maxl,maxm) real pc(maxl,maxm) real psi(maxl,maxm) real rueta(maxl,maxm) real ruksi(maxl,maxm) real t(maxl,maxm) real te(maxl,maxm) real tk(maxl,maxm) real u(maxl,maxm) real v(maxl,maxm) real vmag(maxl,maxm) real vort(maxl,maxm) logical vpflag(maxl,maxm) real w(maxl,maxm) real xbot(maxbot) real xc(maxl,maxm) real xp(maxl,maxm) real ybot(maxbot) real yc(maxl,maxm) real yp(maxl,maxm) ! if ( 2*m-1 > maxm ) then write ( *, * ) ' ' write ( *, * ) 'DOUBLE - Fatal error!' write ( *, * ) ' The reflected data has M = ',2*m-1 write ( *, * ) ' but CPLOT has MAXM = ',maxm stop end if call rduble(e,l,m,maxl,maxm,'e') call rduble(gamt,l,m,maxl,maxm,'gamt') call lduble(ncflag,l,m,maxl,maxm) call lduble(npflag,l,m,maxl,maxm) call rduble(p,l,m,maxl,maxm,'p') call rduble(pc,l,m,maxl,maxm,'pc') call rduble(psi,l,m,maxl,maxm,'psi') call rduble(rueta,l,m,maxl,maxm,'rueta') call rduble(ruksi,l,m,maxl,maxm,'ruksi') call rduble(t,l,m,maxl,maxm,'t') call rduble(te,l,m,maxl,maxm,'te') call rduble(tk,l,m,maxl,maxm,'tk') call rduble(v,l,m,maxl,maxm,'v') call rduble(u,l,m,maxl,maxm,'u') call rduble(vort,l,m,maxl,maxm,'vort') call lduble(vpflag,l,m,maxl,maxm) call rduble(w,l,m,maxl,maxm,'w') call rduble(yp,l,m,maxl,maxm,'yp') call rduble(yc,l,m,maxl,maxm,'yc') call rduble(xp,l,m,maxl,maxm,'xp') call rduble(xc,l,m,maxl,maxm,'xc') call iduble(kcrys,l,m,maxl,maxm) call iduble(kmelt,l,m,maxl,maxm) call iduble(kvoid,l,m,maxl,maxm) ! ! Take care of data defining the crystal extent. ! j = jcrys2 jcrys1 = m-j+1 jcrys2 = m+j-1 ! ! Update M to its new value. ! m = 2*m-1 ! ! Take care of data associated with the crucible shape. ! do i = 2*nbot-1,nbot,-1 xbot(i) = xbot(i+1-nbot) ybot(i) = ybot(i+1-nbot) end do do i = 1,nbot-1 xbot(i) = -xbot(2*nbot-i) ybot(i) = ybot(2*nbot-i) end do nbot = 2*nbot-1 ! ! Set the velocity magnitudes. ! do i = 1,l do j = 1,m vmag(i,j) = sqrt(u(i,j)**2+v(i,j)**2) end do end do return end subroutine dshlin(n,x,y,dshsiz) ! !*****************************************************************************80 ! !! DSHLIN draws a dashed line connecting a series of points. ! ! ! Discussion: ! ! If the X and Y coordinates use different scale ! factors, then dashes at different angles will seem to ! have different lengths. ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of points to be connected. ! ! Input, real X(N), Y(N), the X and Y coordinates of the ! points. ! ! Input, real DSHSIZ, the length, in the X, Y coordinate ! system, of the dashed lines. If it is negative or zero, ! an error occurs. ! implicit none ! integer n ! real dist real dshsiz integer i integer j integer ndash real x(n) real xnew real xold real xxnew real xxold real y(n) real ynew real yold real yynew real yyold ! ! Make sure that DSHSIZ is positive. ! if ( dshsiz <= 0.0 ) then write ( *, * ) ' ' write ( *, * ) 'DshLin - Fatal error!' write ( *, * ) ' The parameter DSHSIZ must be positive.' write ( *, * ) ' but the input value is ',dshsiz return end if ! xnew = x(1) ynew = y(1) do i = 2,n xold = xnew yold = ynew xnew = x(i) ynew = y(i) dist = sqrt( (xnew-xold)**2 + (ynew-yold)**2 ) if ( dist > 0.0 ) then ndash = int((dist/dshsiz)+1) if ( mod(ndash,4) /= 0 ) then ndash = ndash+(4-mod(ndash,4)) end if if ( ndash <= 3)ndash = 4 ! ! Desired pattern is: ! ! X0 - dash - blank - blank - dash - dash - blank - blank - dash - X1 ! do j = 1,ndash if ( mod(j,4) == 0 .or. mod(j,4) == 1 ) then xxold = ( (ndash+1-j)*xold + (j-1)*xnew) / real(ndash) yyold = ( (ndash+1-j)*yold + (j-1)*ynew) / real(ndash) xxnew = ( (ndash-j)*xold + j*xnew) / real(ndash) yynew = ( (ndash-j)*yold + j*ynew) / real(ndash) call movcgm(xxold,yyold) call drwcgm(xxnew,yynew) end if end do end if end do return end subroutine flushl(string) ! !*****************************************************************************80 ! !! FLUSHL flushes a string left. ! ! For instance: ! ! Input Output ! ! ' Hello' 'Hello ' ! ' Hi there! ' 'Hi there! ' ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! STRING Input/output, character*(*) STRING. ! ! On input, STRING is a string of characters. ! ! On output, any initial blank characters in STRING ! have been cut, and pasted back onto the end. ! implicit none ! integer i integer lchar integer nonb character null character ( len = * ) string ! ! Check the length of the string to the last nonblank. ! If nonpositive, return. ! lchar = len_trim (string) if ( lchar <= 0)return null = char(0) ! ! Find the occurrence of the first nonblank. ! do i = 1,lchar nonb = i if ( string(i:i) /= ' ' .and. string(i:i) /= null)go to 10 end do return 10 continue ! ! Shift the string left. ! do i = 1,lchar+1-nonb string(i:i) = string(i+nonb-1:i+nonb-1) end do ! ! Blank out the end of the string. ! do i = lchar+2-nonb,lchar string(i:i) = ' ' end do return end subroutine fsize(l,m,maxl,maxm,nflag,r,rmax,rmin) ! !*****************************************************************************80 ! !! FSIZE computes the range of a real array defined at primary nodes. ! ! ! Only "flagged" values are considered. ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer L, M, the number of rows and columns of data. ! ! Input, integer MAXL, MAXM, the maximum allowed number of ! rows and columns of data. ! ! Input, logical NFLAG(MAXL,MAXM). ! NFLAG is used to "flag" which nodes are active, ! that is, to be displayed, and which not, in the graph. ! ! Input, real R(MAXL,MAXM), the array to be checked. ! ! Output, real RMAX, RMIN, the maximum and minimum values ! of R. ! implicit none ! integer maxl integer maxm ! integer i integer j integer l integer m logical nflag(maxl,maxm) real r(maxl,maxm) real rmax real rmin logical started ! started = .false. rmax = 0.0 rmin = 0.0 do i = 1,l do j = 1,m if ( nflag(i,j) ) then if ( started ) then rmin = min(rmin,r(i,j)) rmax = max(rmax,r(i,j)) else started = .true. rmin = r(i,j) rmax = r(i,j) end if end if end do end do return end subroutine gettab(dev,echo,filgrf,grace,icmax,icmin,ierror,iplot,itable, & lgopen,ovrlay) ! !*****************************************************************************80 ! !! GETTAB gets the color table choice from the user. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character*10 DEV, the graphics output device to be used. ! Current legal values for DEV include: ! ! cgmb - CGM binary file. ! ps - PostScript file. ! xws - X window screen (interactive). ! ! Input, character*80 FILGRF, the name of the output ! graphics file. ! ! Input, real GRACE. ! The size of the "grace" margin on the plot. ! ! Input, integer ICMAX, ICMIN, the maximum and minimum color ! indices to use in color contour graphics. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IPLOT, the number of plots made so far. ! ! Output, integer ITABLE, the desired color table. ! 1: low black to high white ! 2: low blue to high yellow ! 3: low red, high blue, with bands between. ! 4: low red, yellow, green, blue, high white. ! 5: low white, blue, green, yellow, high red. ! 6: low blue to high red ! 7: linear table between 2 user colors. ! 8: linear table between N user colors. ! 9: low white to high black. ! ! Input, logical OVRLAY. ! If OVRLAY is true, then the next time that a plot is ! requested, a "new frame" command is suppressed, so that ! the new plot is shown on top of the previous one. ! implicit none ! character ( len = 10 ) dev logical echo character ( len = 80 ) filcol character ( len = 80 ) filgrf real grace integer icmax integer icmin integer ierror integer iplot integer itable logical lgopen logical ovrlay real x1max real x1min real y1max real y1min ! if ( itable == 0 ) then write ( *, * ) 'Enter name of color table file.' read(*,'(a)',err = 90,end=90)filcol write(17,'(a)')filcol if ( echo ) then write(*,'(a)')filcol end if call getctb(icmin,icmax,filcol,ierror) end if call preplt(dev,echo,filgrf,icmax,icmin,iplot,itable,lgopen,ovrlay) call settab(echo,icmax,icmin,itable) if ( dev == 'xws' ) then call cbox(grace,x1max,x1min,y1max,y1min) call buzz(dev,x1max,x1min,y1max,y1min) end if return 90 continue write ( *, * ) ' ' write ( *, * ) 'GetTab - Error' ierror = 1 return end subroutine getwin(echo,grace,srange,xmax,xmin,x1max,x1min,x2max, & x2min,xsmax,xsmin,ymax,ymin,y1max,y1min,y2max,y2min,ysmax,ysmin) ! !*****************************************************************************80 ! !! GETWIN responds to the "W" command by telling the user ! the current window, and getting the new value from the ! user. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real GRACE. ! The size of the "grace" margin on the plot. ! ! SRANGE Output, real SRANGE. ! The maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! XMAX Input, real XMAX. ! The maximum X coordinate of all the nodes. ! The maximum entry in the XC array. ! ! XMIN Input, real XMIN. ! The minimum X coordinate of all the nodes. ! The minimum entry in the XC array. ! ! X1MAX, ! X1MIN Output, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! X2MAX, ! X2MIN Input/output, real X2MAX, X2MIN, the maximum and minimum X ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! XSMAX Input/output, real XSMAX. ! The maximum X coordinate of the data to be displayed. ! XSMAX defaults to XMAX, but can be made smaller to ! focus on a portion of the region. ! ! XSMIN Input/output, real XSMIN. ! The minimum X coordinate of the data to be displayed. ! XSMIN defaults to XMIN, but can be made larger to ! focus on a portion of the region. ! ! Y1MAX, ! Y1MIN Output, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! ! Y2MAX, ! Y2MIN Input/output, real Y2MAX, Y2MIN, the maximum and minimum Y ! coordinates that should be used for plotting. No plotting commands ! should exceed these values. This is where the "frame" might be ! drawn. ! ! YMAX Input, real YMAX. ! The maximum Y coordinate of all the nodes. ! The maximum value attained by the YC array. ! ! YMIN Input, real YMIN. ! The minimum Y coordinate of all the nodes. ! The minimum value attained by the YC array. ! ! YSMAX Input/output, real YSMAX. ! The maximum Y coordinate of the data to be displayed. ! YSMAX defaults to YMAX, but can be made smaller to ! focus on a portion of the region. ! ! YSMIN Input/output, real YSMIN. ! The minimum Y coordinate of the data to be displayed. ! YSMIN defaults to YMIN, but can be made larger to ! focus on a portion of the region. ! implicit none ! logical echo real grace real srange real x1max real x1min real x2max real x2min real xmax real xmin real xsmax real xsmin real y1max real y1min real y2max real y2min real ymax real ymin real ysmax real ysmin ! 10 continue write ( *, * ) ' ' write ( *, * ) 'GetWin:' write ( *, * ) ' ' write ( *, * ) ' Data coordinates:' write ( *, * ) ' ' write ( *, * ) xmin,' = XMIN <= X <= XMAX = ',xmax write ( *, * ) ymin,' = YMIN <= Y <= YMAX = ',ymax write ( *, * ) ' ' write ( *, * ) ' Displayed data coordinates:' write ( *, * ) ' ' write ( *, * ) xsmin,' = XSMIN <= X <= XSMAX = ',xsmax write ( *, * ) ysmin,' = YSMIN <= Y <= YSMAX = ',ysmax write ( *, * ) ' ' write ( *, * ) ' Picture window coordinates:' write ( *, * ) ' ' write ( *, * ) x2min,' = X2MIN <= X <= X2MAX = ',x2max write ( *, * ) y2min,' = Y2MIN <= Y <= Y2MAX = ',y2max write ( *, * ) ' ' write ( *, * ) ' Enter new displayed data coordinates for X and Y:' write ( *, * ) ' Use the order XSMIN, XSMAX, YSMIN, YSMAX:' read(*,*,err = 20,end=20)xsmin,xsmax,ysmin,ysmax write(17,*)xsmin,xsmax,ysmin,ysmax if ( echo ) then write ( *, * ) xsmin,xsmax,ysmin,ysmax end if ! ! Compute box containing data. ! call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) return 20 continue write ( *, * ) ' ' write ( *, * ) 'GetWin - Warning!' write ( *, * ) ' Your input was not acceptable!' write ( *, * ) ' Please try again!' write ( *, * ) ' ' go to 10 end subroutine graph(delx,dely,dev,dxcdp,dycdp,echo,filgrf,gamt,icmax,icmin, & icolor,icrys1,icrys2,iplot,itable,jcmax,jcmin,jcrys1,jcrys2,kcrys,kmelt, & kvoid,l,lbar,lflag,lgopen,m,maxbot,maxl,maxm,maxobj,nbot,ncflag,ncon, & nflag,npflag,nxskip,nyskip,object,ovrlay,p,pc,psi,s1,s2,scalecv,scalenb, & scalenc,scalenp,scalev,show,srange,t,te,title,title2,tk,u,v,vpflag,vmag, & vort,w,x1max,x1min,x2max,x2min,xbot,xc,xp,xsmax,xsmin,y1max,y1min, & y2max,y2min,ybot,yc,yp,ysmax,ysmin) ! !*****************************************************************************80 ! !! GRAPH draws the graph, after the user has specified what is to be shown. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real DELX. ! The X spacing between nodes. In some cases, ! this spacing is modified to create isoparametric elements. ! ! DELY Input, real DELY. ! The Y spacing between nodes. In some cases, ! this spacing is modified to create isoparametric elements. ! ! DEV Input, character*10 DEV. ! The graphics output device to be used. Current legal ! values for DEV include: ! ! cgmb - CGM binary file. ! ps - PostScript file. ! xws - X window screen (interactive). ! ! DXCDP Input, real DXCDP(MAXL,MAXM), the difference between ! the second and first values of XC. ! ! DYCDP Input, real DYCDP(MAXL,MAXM), the difference between ! the second and first values of YC. ! ! FILGRF Input, character*80 FILGRF, the name of the output ! graphics file. ! ! GAMT Output, real GAMT(MAXL,MAXM), the diffusion coefficient. ! ! ICMAX, ! ICMIN Input, integer ICMAX, ICMIN, the maximum and ! minimum color indices to use in the color bar. ! ! ICOLOR Input, integer ICOLOR(MAXOBJ). ! Contains the color indexes for each object. ! However, in some cases, ICOLOR is actual a color table ! index. ! ! ICRYS1 Input, integer ICRYS1, the I coordinate of the row of corner ! nodes that define the bottom of the crystal. ! ! ICRYS2 Input, integer ICRYS2, the I coordinate of the row of corner ! nodes that define the top of the crystal. ! ! IPLOT Input, integer IPLOT. ! The number of plots made so far. ! ! ITABLE Input, integer ITABLE, the desired color table. ! ! 1: low black to high white ! 2: low blue to high yellow ! 3: low red, high blue, with bands between. ! 4: low red, yellow, green, blue, high white. ! 5: low white, blue, green, yellow, high red. ! 6: low blue to high red ! 7: linear table between 2 user colors. ! 8: linear table between N user colors. ! 9: low white to high black. ! ! ! JCMAX, ! JCMIN Input, integer JCMAX, JCMIN, the maximum and ! minimum color indices to use in the color bar. ! ! JCRYS1 Input/output, integer JCRYS1, the J coordinate of the column of ! corner nodes that define the left side of the crystal. ! ! JCRYS2 Input/output, integer JCRYS2, the J coordinate of the column of ! corner nodes that define the right side of the crystal. ! ! KCRYS Input, integer KCRYS(MAXL,MAXM), ! 0, if control volume (I,J) is away from the crystal. ! 1, if control volume (I,J) is on the external crystal boundary. ! 2, if control volume (I,J) is on the internal crystal boundary. ! 3, if control volume (I,J) is in the crystal interior. ! ! LFLAG Workspace, logical LFLAG(MAXL,MAXM), .TRUE. if primary node ! (I,J) lies within the visible data, .FALSE. otherwise. ! ! KMELT Input, integer KMELT(MAXL,MAXM), ! 0, if control volume (I,J) is away from the melt. ! 1, if control volume (I,J) is on the external melt boundary. ! 2, if control volume (I,J) is on the internal melt boundary. ! 3, if control volume (I,J) is in the melt interior. ! ! KVOID Input, integer KVOID(MAXL,MAXM), ! 0, if control volume (I,J) is away from the void. ! 1, if control volume (I,J) is on the external void boundary. ! 2, if control volume (I,J) is on the internal void boundary. ! 3, if control volume (I,J) is in the void interior. ! ! L Input, integer L, the number of rows of data. ! ! LBAR Input, logical LBAR, is .TRUE. if the color bar should ! be shown. ! ! M Input, integer M, the number of columns of data. ! ! MAXBOT Input, integer MAXBOT, the maximum number of crucible nodes ! allowed. ! ! MAXL, ! MAXM Input, integer MAXL, MAXM, the maximum allowed number of ! rows and columns of data. ! ! MAXOBJ Input, integer MAXOBJ. ! The number of graphical "objects". ! ! NBOT Input, integer NBOT, the number of crucible nodes. ! ! NCON Input, integer NCON. ! The number of contour lines to be drawn. This is ! initialized to 12, but may be changed by the user. ! ! NFLAG Input, logical NFLAG(MAXL,MAXM). ! ! NFLAG is used to "flag" which nodes are active, ! that is, to be displayed, and which not, in the graph. ! ! NXSKIP Input, integer NXSKIP. ! NXSKIP is used to "thin" out a vector plot. ! ! If NXSKIP = 1, then a standard vector plot is made. ! ! Otherwise, in the X direction, vectors are drawn only ! in columns 1, 1+NXSKIP, 1+2*NXSKIP and so on. ! ! NYSKIP Input, integer NYSKIP. ! NYSKIP is used to "thin" out a vector plot. ! ! If NYSKIP = 1, then a standard vector plot is made. ! ! Otherwise, in the Y direction, vectors are drawn only ! in rows 1, 1+NYSKIP, 1+2*NYSKIP and so on. ! ! OBJECT Input, character*25 OBJECT(MAXOBJ), the names of the ! graphical objects. ! ! OVRLAY Input, logical OVRLAY. ! If OVRLAY is true, then the next time that a plot is ! requested, a "new frame" command is suppressed, so that ! the new plot is shown on top of the previous one. ! ! P Input, real P(MAXL,MAXM), the pressure. ! ! PC Input, real PC(MAXL,MAXM), the corrected pressure. ! ! PSI Input, real PSI(MAXL,MAXM), the stream function. ! ! S1 Workspace, real S1(MAXL,MAXM). ! ! S2 Workspace, real S2(MAXL,MAXM). ! ! SHOW Input, logical SHOW(MAXOBJ). ! Contains, for each object, a flag determining whether it ! is to be shown or not. ! ! SRANGE Input, real SRANGE. ! The maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! SCALECV Input, real SCALECV. ! A scale factor for the control volume numbers, which has a ! default value of 1. ! ! SCALENB Input, real SCALENB. ! A scale factor for the size of the nodes that define the ! crucible shape, which has a default value of 1. ! ! SCALENC Input, real SCALENC. ! A scale factor for the size of the corner nodes, which has a ! default value of 1. ! ! SCALENP Input, real SCALENP. ! A scale factor for the size of the primary nodes, which has a ! default value of 1. ! ! SCALEV Input, real SCALEV. ! A scale factor for velocity vectors. This starts out at 1.0. ! ! T Input, real T(MAXL,MAXM), the temperature. ! ! TE Input, real TE(MAXL,MAXM), the turbulent epsilon. ! ! TITLE Input, character*40 TITLE. ! A title for the plots. ! ! TITLE2 Input, character*40 TITLE2. ! A subtitle used in the profile plots. ! ! TK Input, real TK(MAXL,MAXM), the turbulent K. ! ! U Input, real U(MAXL,MAXM), the horizontal velocity. ! ! V Input, real V(MAXL,MAXM), the vertical velocity. ! ! VPFLAG Input, logical VPFLAG(MAXL,MAXM), a flag which determines ! whether vector data at a given node should be displayed. ! ! VMAG Input, real VMAG(MAXL,MAXM), the velocity magnitude. ! ! W Input, real W(MAXL,MAXM), the axial velocity. ! ! X Output, real X(MAXL,MAXM), the X coordinates of the primary ! nodes. ! ! X1MAX, ! X1MIN Output, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! X2MAX, ! X2MIN Input, real X2MAX, X2MIN, the maximum and minimum X ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! XBOT Input, real XBOT(NBOT), the X coordinates of the crucible ! bottom nodes. ! ! XC Input, real XC(MAXL,MAXM), the X coordinates of the corner ! nodes. ! ! XSMAX Input, real XSMAX. ! The maximum X coordinate of the data to be displayed. ! XSMAX defaults to XMAX, but can be made smaller to ! focus on a portion of the region. ! ! XSMIN Input, real XSMIN. ! The minimum X coordinate of the data to be displayed. ! XSMIN defaults to XMIN, but can be made larger to ! focus on a portion of the region. ! ! Y Output, real Y(MAXL,MAXM), the Y coordinates of the primary ! nodes. ! ! Y1MAX, ! Y1MIN Output, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! ! Y2MAX, ! Y2MIN Input, real Y2MAX, Y2MIN, the maximum and minimum Y ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! YBOT Input, real YBOT(NBOT), the Y coordinates of the crucible ! bottom nodes. ! ! YC Input, real YC(MAXL,MAXM), the Y coordinates of the corner ! nodes. ! ! YSMAX Input, real YSMAX. ! The maximum Y coordinate of the data to be displayed. ! YSMAX defaults to YMAX, but can be made smaller to ! focus on a portion of the region. ! ! YSMIN Input, real YSMIN. ! The minimum Y coordinate of the data to be displayed. ! YSMIN defaults to YMIN, but can be made larger to ! focus on a portion of the region. ! implicit none ! integer maxbot integer maxl integer maxm integer maxobj ! character ( len = 25 ) aname real angle character ( len = 6 ) chrtmp real csize real cwide real delx real dely character ( len = 10 ) dev real dshsiz real dxcdp(maxl,maxm) real dycdp(maxl,maxm) logical echo character ( len = 80 ) filgrf logical filled character flush real gamt(maxl,maxm) integer i integer icmax integer icmin integer icolor(maxobj) integer icrys1 integer icrys2 integer ido logical inside integer iplot integer itable integer j integer jcmax integer jcmin integer jcolor integer jcrys1 integer jcrys2 integer kcrys(maxl,maxm) integer kmelt(maxl,maxm) integer kvoid(maxl,maxm) integer l logical lbar logical lbar2 integer lent logical lflag(maxl,maxm) logical lgopen integer m integer nbot logical ncflag(maxl,maxm) integer ncon logical nflag(maxl,maxm) integer nonzer logical npflag(maxl,maxm) integer npts integer nval integer nxskip integer nyskip character ( len = 25 ) object(maxobj) logical ovrlay real p(maxl,maxm) real pc(maxl,maxm) real psi(maxl,maxm) real pwide real s1(maxl,maxm) real s2(maxl,maxm) real scalecv real scalenb real scalenc real scalenp real scalev logical show(maxobj) real srange real t(maxl,maxm) real te(maxl,maxm) character ( len = 40 ) title character ( len = 40 ) title2 real tk(maxl,maxm) real u(maxl,maxm) real uu real uvmag real v(maxl,maxm) real velscl logical vpflag(maxl,maxm) real vmag(maxl,maxm) real vmax real vort(maxl,maxm) real vv real w(maxl,maxm) real x1max real x1min real x2max real x2min real xbot(maxbot) real xc(maxl,maxm) real xp(maxl,maxm) real xpoly(4) real xsmax real xsmin real xval real xvec(2) real xx(300) real y1max real y1min real y2max real y2min real ybot(maxbot) real yc(maxl,maxm) real yp(maxl,maxm) real ypoly(4) real ysmax real ysmin real yval real yvec(2) real yy(300) ! call preplt(dev,echo,filgrf,icmax,icmin,iplot,itable,lgopen,ovrlay) ! cwide = 0.9*srange/40 pwide = srange ! ! Set the scale of the picture. ! We allow ourselves more or less a one percent margin. ! xvec(1) = x2min-0.01*(x2max-x2min) yvec(1) = y2min-0.01*(y2max-y2min) xvec(2) = x2max+0.01*(x2max-x2min) yvec(2) = y2max+0.01*(y2max-y2min) call setscl(xvec,yvec,2) ! ! Mark those primary nodes which are visible. ! ! This marker is based on visibility (inside the screen limits) ! and on whether the user has asked to hide the crystal, ! melt, or void regions. ! do i = 1,l do j = 1,m if ( (.not.inside(xsmin,xp(i,j),xsmax)) .or. & (.not.inside(ysmin,yp(i,j),ysmax)) ) then lflag(i,j) = .false. else if ( kcrys(i,j) == 3 .and. (.not.show(39)) ) then lflag(i,j) = .false. else if ( kmelt(i,j) == 3 .and. (.not.show(40)) ) then lflag(i,j) = .false. else if ( kvoid(i,j) == 3 .and. (.not.show(38)) ) then lflag(i,j) = .false. else lflag(i,j) = .true. end if end do end do ! ! FRAME: Draw a box around the region. ! if ( show(3) ) then call linclr(icolor(3)) call box(x2min,x2max,y2min,y2max) end if ! ! TITLE/TITLE2: Draw the titles. ! if ( show(7) ) then call linclr(icolor(7)) lent = len_trim ( title ) if ( lent > 0 ) then angle = 0.0 xval = 0.5*(xsmax+xsmin) yval = ysmax+3.0*cwide flush = 'center' call s_plot(angle,cwide,pwide,title,xval,yval,flush) end if lent = len_trim ( title2 ) if ( lent > 0 ) then angle = 0.0 xval = 0.5*(xsmax+xsmin) yval = ysmax+1.5*cwide flush = 'center' call s_plot(angle,cwide,pwide,title2,xval,yval,flush) end if end if ! ! CPC: Draw corrected pressure color contours for nodes which are visible, ! and which are in or on the melt region. ! if ( show(26) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do aname = object(26) call colcon(pc,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! CRYSC: color in the crystal. ! if ( show(44) ) then jcolor = icolor(44) call filclr(jcolor) do i = 1,l-1 do j = 1,m-1 if ( kcrys(i,j) >= 2 ) then npts = 4 xpoly(1) = xc(i,j) ypoly(1) = yc(i,j) xpoly(2) = xc(i+1,j) ypoly(2) = yc(i+1,j) xpoly(3) = xc(i+1,j+1) ypoly(3) = yc(i+1,j+1) xpoly(4) = xc(i,j+1) ypoly(4) = yc(i,j+1) call plygon(npts,xpoly,ypoly) end if end do end do end if ! ! MELTC: color in the melt. ! if ( show(45) ) then jcolor = icolor(45) call filclr(jcolor) do i = 1,l-1 do j = 1,m-1 if ( kmelt(i,j) >= 2 ) then npts = 4 xpoly(1) = xc(i,j) ypoly(1) = yc(i,j) xpoly(2) = xc(i+1,j) ypoly(2) = yc(i+1,j) xpoly(3) = xc(i+1,j+1) ypoly(3) = yc(i+1,j+1) xpoly(4) = xc(i,j+1) ypoly(4) = yc(i,j+1) call plygon(npts,xpoly,ypoly) end if end do end do end if ! ! PC: Draw pressure color contours for nodes which are visible, ! and which are in or on the melt region. ! if ( show(12) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do aname = object(12) call colcon(p,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! SC: Draw stream color contours in the melt region. ! if ( show(35) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do aname = object(35) call colcon(psi,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! TEMPC: Draw temperature color contours for nodes which are visible, ! and which are in or on the melt or crystal regions. ! if ( show(13) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. (kmelt(i,j) >= 2 .or. kcrys(i,j) >= 2) ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do aname = object(13) call colcon(t,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! VMAGC: Draw velocity magnitude color contours for nodes which are visible, ! and which are in or on the melt region. ! if ( show(37) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do aname = object(37) call colcon(vmag,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! VOIDC: color in the void. ! if ( show(46) ) then jcolor = icolor(46) call filclr(jcolor) do i = 1,l-1 do j = 1,m-1 if ( kvoid(i,j) >= 2 ) then npts = 4 xpoly(1) = xc(i,j) ypoly(1) = yc(i,j) xpoly(2) = xc(i+1,j) ypoly(2) = yc(i+1,j) xpoly(3) = xc(i+1,j+1) ypoly(3) = yc(i+1,j+1) xpoly(4) = xc(i,j+1) ypoly(4) = yc(i,j+1) call plygon(npts,xpoly,ypoly) end if end do end do end if ! ! VORTC: draw vorticity color contours in melt region. ! if ( show(42) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do aname = object(42) call colcon(vort,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! WC: Draw angular momentum color contours in melt region. ! if ( show(14) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do aname = object(14) call colcon(w,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! XCC: draw XC coordinate color contours. ! if ( show(31) ) then do i = 1,l do j = 1,m nflag(i,j) = ncflag(i,j) if ( .not.lflag(i,j) ) then nflag(i,j) = .false. end if end do end do aname = object(31) call colcon(xc,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! YCC: draw YC coordinate color contours. ! if ( show(33) ) then do i = 1,l do j = 1,m nflag(i,j) = ncflag(i,j) if ( .not.lflag(i,j) ) then nflag(i,j) = .false. end if end do end do aname = object(33) call colcon(yc,aname,echo,icolor,jcmax,jcmin,l,lbar,m,maxl, & maxm,maxobj,ncon,nflag,srange,xp,x2max,x2min,yp,y1max,y2max) end if ! ! LINE DRAWINGS BEGIN HERE. ! ! ! AX: draw the axis of symmetry. ! if ( show(28) ) then if ( xsmin <= 0.0 .and. 0.0 <= xsmax ) then call linclr(icolor(28)) nval = 2 xx(1) = 0.0 yy(1) = ysmin xx(2) = 0.0 yy(2) = ysmax dshsiz = 0.025 call dshlin(nval,xx,yy,dshsiz) end if end if ! ! BO: Draw the boundary. ! if ( show(1) ) then call linclr(icolor(1)) dshsiz = 0.05 nval = 2 do i = 1,l-1 if ( ncflag(i,1) .and. ncflag(i+1,1) ) then xx(1) = xc(i,1) yy(1) = yc(i,1) xx(2) = xc(i+1,1) yy(2) = yc(i+1,1) call dshlin(nval,xx,yy,dshsiz) end if end do do j = 1,m-1 if ( ncflag(l,j) .and. ncflag(l,j+1) ) then xx(1) = xc(l,j) yy(1) = yc(l,j) xx(2) = xc(l,j+1) yy(2) = yc(l,j+1) call dshlin(nval,xx,yy,dshsiz) end if end do do i = 1,l-1 if ( ncflag(i,m) .and. ncflag(i+1,m) ) then xx(1) = xc(i,m) yy(1) = yc(i,m) xx(2) = xc(i+1,m) yy(2) = yc(i+1,m) call dshlin(nval,xx,yy,dshsiz) end if end do do j = 1,m-1 if ( ncflag(1,j) .and. ncflag(1,j+1) ) then xx(1) = xc(1,j) yy(1) = yc(1,j) xx(2) = xc(1,j+1) yy(2) = yc(1,j+1) call dshlin(nval,xx,yy,dshsiz) end if end do end if ! ! CP: Draw corrected pressure contours in melt region. ! if ( show(25) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(26) ) then ido = 0 call linclr(icolor(25)) else ido = 1 end if aname = object(25) lbar2 = lbar .and. .not.show(26) call lincon(pc,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! CRU: Draw the crucible boundary. ! if ( show(19) ) then call linclr(icolor(19)) i = 1 do j = 1,m-1 if ( xsmin <= xc(i,j) .and. xc(i,j) <= xsmax.and. & ysmin <= yc(i,j) .and. yc(i,j) <= ysmax.and. & xsmin <= xc(i,j+1) .and. xc(i,j+1) <= xsmax.and. & ysmin <= yc(i,j+1) .and. yc(i,j+1) <= ysmax ) then call movcgm(xc(i,j),yc(i,j)) call drwcgm(xc(i,j+1),yc(i,j+1)) end if end do ! ! Do this loop only if reflected. ! if ( xc(1,1) == -xc(1,m) ) then j = 1 do i = 1,icrys1-1 if ( xsmin <= xc(i,j) .and. xc(i,j) <= xsmax.and. & ysmin <= yc(i,j) .and. yc(i,j) <= ysmax.and. & xsmin <= xc(i+1,j) .and. xc(i+1,j) <= xsmax.and. & ysmin <= yc(i+1,j) .and. yc(i+1,j) <= ysmax ) then call movcgm(xc(i,j),yc(i,j)) call drwcgm(xc(i+1,j),yc(i+1,j)) end if end do end if j = m do i = 1,icrys1-1 if ( xsmin <= xc(i,j) .and. xc(i,j) <= xsmax.and. & ysmin <= yc(i,j) .and. yc(i,j) <= ysmax.and. & xsmin <= xc(i+1,j) .and. xc(i+1,j) <= xsmax.and. & ysmin <= yc(i+1,j) .and. yc(i+1,j) <= ysmax ) then call movcgm(xc(i,j),yc(i,j)) call drwcgm(xc(i+1,j),yc(i+1,j)) end if end do end if ! ! CRYS: Draw the crystal. ! if ( show(20) ) then call linclr(icolor(20)) call movcgm(xc(icrys1,jcrys1),yc(icrys1,jcrys1)) do j = jcrys1+1,jcrys2 if ( ncflag(icrys1,j-1) .and. ncflag(icrys1,j) ) then call drwcgm(xc(icrys1,j),yc(icrys1,j)) end if end do do i = icrys1+1,icrys2 if ( ncflag(i-1,jcrys2) .and. ncflag(i,jcrys2) ) then call drwcgm(xc(i,jcrys2),yc(i,jcrys2)) end if end do do j = jcrys2-1,jcrys1,-1 if ( ncflag(icrys2,j+1) .and. ncflag(icrys2,j) ) then call drwcgm(xc(icrys2,j),yc(icrys2,j)) end if end do do i = icrys2,icrys1,-1 if ( ncflag(i+1,jcrys1) .and. ncflag(i,jcrys1) ) then call drwcgm(xc(i,jcrys1),yc(i,jcrys1)) end if end do end if ! ! CV: Draw the control volumes. ! if ( show(2) ) then call linclr(icolor(2)) do i = 1,l do j = 1,m if ( (.not.inside(xsmin,xc(i,j),xsmax)) .or. & (.not.inside(ysmin,yc(i,j),ysmax)) ) then nflag(i,j) = .false. else nflag(i,j) = .true. end if end do end do do i = 1,l-1 do j = 1,m-1 if ( (kcrys(i,j) == 3 .and. (.not.show(39))) .or. & (kmelt(i,j) == 3 .and. (.not.show(40))) .or. & (kvoid(i,j) == 3 .and. (.not.show(38))) ) then nflag(i,j) = .false. nflag(i+1,j) = .false. nflag(i,j+1) = .false. nflag(i+1,j+1) = .false. end if end do end do ! do i = 1,l ! do j = 1,m ! nflag(i,j) = lflag(i,j) ! end do ! end do call showcv(l,m,maxl,maxm,nflag,xc,yc) end if ! ! CVN: Draw the control volume numbers. ! if ( show(43) ) then angle = 0.0 cwide = 0.05*srange*scalecv pwide = srange flush = 'center' call linclr(icolor(43)) do i = 1,l-1 do j = 1,m-1 if ( lflag(i,j) ) then write(chrtmp,'(i6)')i call flushl(chrtmp) call s_plot(angle,cwide,pwide,chrtmp,xp(i,j),yp(i,j),flush) end if end do end do end if ! ! FS: Draw the free surface. ! if ( show(21) ) then call linclr(icolor(21)) do i = 2,l-1 do j = 1,m-1 if ( kmelt(i,j) == 1 .and. kmelt(i-1,j) == 2) then if ( xsmin <= xc(i,j) .and. xc(i,j) <= xsmax .and. & ysmin <= yc(i,j) .and. yc(i,j) <= ysmax .and. & xsmin <= xc(i,j+1) .and. xc(i,j+1) <= xsmax .and. & ysmin <= yc(i,j+1) .and. yc(i,j+1) <= ysmax ) then call movcgm(xc(i,j),yc(i,j)) call drwcgm(xc(i,j+1),yc(i,j+1)) end if end if end do end do end if ! ! NB: Draw the crucible nodes. ! if ( show(34) ) then call linclr(icolor(34)) filled = .false. csize = 0.0125*srange*scalenb do i = 1,nbot if ( inside(xsmin,xbot(i),xsmax) .and. & inside(ysmin,ybot(i),ysmax) ) then call circle(xbot(i),ybot(i),csize,filled) end if end do end if ! ! NC: Draw the corner nodes. ! if ( show(29) ) then call linclr(icolor(29)) filled = .false. csize = 0.0125*srange*scalenc do i = 1,l do j = 1,m nflag(i,j) = ncflag(i,j) if ( mod(i-1,nxskip) /= 0 ) then nflag(i,j) = .false. end if if ( mod(j-1,nyskip) /= 0 ) then nflag(i,j) = .false. end if if ( nflag(i,j) ) then call diamnd(xc(i,j),yc(i,j),csize,filled) end if end do end do end if ! ! NP: Draw the primary nodes. ! if ( show(4) ) then call linclr(icolor(4)) filled = .false. csize = 0.0125*srange*scalenp do i = 1,l do j = 1,m nflag(i,j) = npflag(i,j) if ( mod(i-1,nxskip) /= 0 ) then nflag(i,j) = .false. end if if ( mod(j-1,nyskip) /= 0 ) then nflag(i,j) = .false. end if if ( nflag(i,j) ) then call circle(xp(i,j),yp(i,j),csize,filled) end if end do end do end if ! ! P: Draw pressure contours in melt region. ! if ( show(5) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(12) ) then ido = 0 call linclr(icolor(5)) else ido = 1 end if aname = object(5) lbar2 = lbar .and. .not.show(12) call lincon(p,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! S: Draw stream lines in melt region. ! if ( show(6) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(35) ) then ido = 0 call linclr(icolor(6)) else ido = 1 end if ! ! Temporary ! ido = 0 aname = object(6) lbar2 = lbar call lincon(psi,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! TE: Draw turbulent epsilon contours in melt region. ! if ( show(22) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do ido = 1 aname = object(22) lbar2 = lbar call lincon(te,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! TEMP: Draw temperature contours in melt and crystal regions. ! if ( show(11) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. (kmelt(i,j) >= 2 .or. kcrys(i,j) >= 2) ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(13) ) then ido = 0 call linclr(icolor(11)) else ido = 1 end if aname = object(11) lbar2 = lbar .and. .not.show(13) call lincon(t,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! TK: Draw turbulent KE contours in melt region. ! if ( show(23) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do ido = 1 aname = object(23) lbar2 = lbar call lincon(tk,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! VIS: Draw viscosity contours in melt region. ! if ( show(24) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do ido = 1 aname = object(24) lbar2 = lbar call lincon(gamt,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! VMAG: Draw velocity magnitude contours in melt region. ! if ( show(36) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(37) ) then ido = 0 call linclr(icolor(36)) else ido = 1 end if aname = object(36) lbar2 = lbar .and. .not.show(37) call lincon(vmag,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! VORT: Draw vorticity contours in melt. ! if ( show(41) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(42) ) then ido = 0 call linclr(icolor(41)) else ido = 1 end if aname = object(41) lbar2 = lbar .and. .not.show(42) call lincon(vort,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! W: Draw angular momentum contours in melt region. ! if ( show(10) ) then do i = 1,l do j = 1,m if ( lflag(i,j) .and. kmelt(i,j) >= 2 ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(14) ) then ido = 0 call linclr(icolor(10)) else ido = 1 end if aname = object(10) lbar2 = lbar .and. .not.show(14) call lincon(w,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! XC: Draw X coordinate contours. ! if ( show(30) ) then do i = 1,l do j = 1,m if ( lflag(i,j) ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(31) ) then ido = 0 call linclr(icolor(30)) else ido = 1 end if aname = object(30) lbar2 = lbar .and. .not.show(31) call lincon(xp,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! YC: Draw Y coordinate contours. ! if ( show(32) ) then do i = 1,l do j = 1,m if ( lflag(i,j) ) then nflag(i,j) = .true. else nflag(i,j) = .false. end if end do end do if ( show(33) ) then ido = 0 call linclr(icolor(32)) else ido = 1 end if aname = object(32) lbar2 = lbar .and. .not.show(33) call lincon(yp,aname,echo,icolor,ido,jcmax,jcmin,l,lbar2,m, & maxl,maxm,maxobj,ncon,nflag,srange,x2max,x2min,xp,y1max,y2max,yp) end if ! ! V: Draw velocity vectors in melt region. ! if ( show(8) ) then call linclr(icolor(8)) do i = 1,l do j = 1,m nflag(i,j) = npflag(i,j) if ( mod(i-1,nxskip) /= 0 ) then nflag(i,j) = .false. end if if ( mod(j-1,nyskip) /= 0 ) then nflag(i,j) = .false. end if if ( .not.vpflag(i,j) ) then nflag(i,j) = .false. end if if ( kmelt(i,j) <= 1 ) then nflag(i,j) = .false. end if if ( .not.lflag(i,j) ) then nflag(i,j) = .false. end if end do end do ! ! Get the maximum velocity norm of the displayed data. ! call vsize(l,m,maxl,maxm,nflag,u,v,vmax) if ( vmax == 0.0 ) then write ( *, * ) ' ' write ( *, * ) 'GRAPH - Warning!' write ( *, * ) ' No nonzero velocity vectors were visible.' else velscl = scalev*0.5*min(nxskip*delx,nyskip*dely)/vmax call vector(l,m,maxl,maxm,nflag,velscl,u,v,xp,yp) end if end if ! ! U: Draw velocity direction vectors in melt region. ! if ( show(9) ) then call linclr(icolor(9)) do i = 1,l do j = 1,m nflag(i,j) = npflag(i,j) if ( mod(i-1,nxskip) /= 0 ) then nflag(i,j) = .false. end if if ( mod(j-1,nyskip) /= 0 ) then nflag(i,j) = .false. end if if ( .not.vpflag(i,j) ) then nflag(i,j) = .false. end if if ( kmelt(i,j) <= 1 ) then nflag(i,j) = .false. end if if ( .not.lflag(i,j) ) then nflag(i,j) = .false. end if end do end do nonzer = 0 do i = 1,l do j = 1,m uu = u(i,j) vv = v(i,j) uvmag = sqrt(uu*uu+vv*vv) if ( uvmag /= 0.0 ) then nonzer = nonzer+1 s1(i,j) = uu/uvmag s2(i,j) = vv/uvmag else s1(i,j) = 0.0 s2(i,j) = 0.0 end if end do end do if ( nonzer > 0 ) then velscl = scalev*0.5*min(nyskip*delx,nxskip*dely) call vector(l,m,maxl,maxm,nflag,velscl,s1,s2,xp,yp) else write ( *, * ) ' ' write ( *, * ) 'GRAPH - Warning!' write ( *, * ) ' No nonzero velocity vectors were visible.' end if end if ! ! Draw (dXCdP, dYCdP) direction vectors in entire region. ! if ( show(27) ) then call linclr(icolor(27)) do i = 1,l do j = 1,m nflag(i,j) = ncflag(i,j) if ( mod(i-1,nxskip) /= 0 ) then nflag(i,j) = .false. end if if ( mod(j-1,nyskip) /= 0 ) then nflag(i,j) = .false. end if if ( .not.lflag(i,j) ) then nflag(i,j) = .false. end if end do end do nonzer = 0 do i = 1,l do j = 1,m uu = dxcdp(i,j) vv = dycdp(i,j) uvmag = sqrt(uu*uu+vv*vv) if ( uvmag /= 0.0 ) then nonzer = nonzer+1 s1(i,j) = uu/uvmag s2(i,j) = vv/uvmag else s1(i,j) = 0.0 s2(i,j) = 0.0 end if end do end do if ( nonzer > 0 ) then velscl = scalev*0.5*min(nyskip*delx,nxskip*dely) call vector(l,m,maxl,maxm,nflag,velscl,s1,s2,xp,yp) else write ( *, * ) ' ' write ( *, * ) 'GRAPH - Warning!' write ( *, * ) ' No (dXCdP,dYCdP) directions, all zero.' end if end if ! ! Pause, if we are doing X-Windows. ! call buzz(dev,x1min,x1max,y1min,y1max) return end subroutine half(e,gamt,jcrys1,jcrys2,kcrys,kmelt,kvoid,l,m, & maxbot,maxl,maxm,nbot,ncflag,npflag,p,pc,psi,rueta,ruksi, & t,te,tk,u,v,vmag,vort,vpflag,w,xbot,xc,xp,ybot,yc,yp) ! !*****************************************************************************80 ! !! HALF removes the reflected data. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real E(MAXL,MAXM), the magnetic stream function. ! ! Input/output, real GAMT(MAXL,MAXM), the diffusion coefficient. ! ! Input/output, integer JCRYS1, the J coordinate of the column of ! corner nodes that define the left side of the crystal. ! ! Input/output, integer JCRYS2, the J coordinate of the column of ! corner nodes that define the right side of the crystal. ! ! Input/output, integer KCRYS(MAXL,MAXM), ! 0, if control volume (I,J) is away from the crystal. ! 1, if control volume (I,J) is on the external crystal boundary. ! 2, if control volume (I,J) is on the internal crystal boundary. ! 3, if control volume (I,J) is in the crystal interior. ! ! Input/output, integer KMELT(MAXL,MAXM), ! 0, if control volume (I,J) is away from the melt. ! 1, if control volume (I,J) is on the external melt boundary. ! 2, if control volume (I,J) is on the internal melt boundary. ! 3, if control volume (I,J) is in the melt interior. ! ! Input/output, integer KVOID(MAXL,MAXM), ! 0, if control volume (I,J) is away from the void. ! 1, if control volume (I,J) is on the external void boundary. ! 2, if control volume (I,J) is on the internal void boundary. ! 3, if control volume (I,J) is in the void interior. ! ! Input, integer L, the number of rows of data. ! ! Input/output, integer M, the number of columns of data. ! ! Input, integer MAXL, MAXM, the maximum allowed number of ! rows and columns of data. ! ! Input/output, real P(MAXL,MAXM), the pressure. ! ! Input/output, real PC(MAXL,MAXM), the corrected pressure. ! ! Input/output, real PSI(MAXL,MAXM), the stream function. ! ! Input/output, real RUETA(MAXL,MAXM), RUKSI(MAXL,MAXM), the ! the momentum in the ETA and KSI directions. ! ! Input/output, real T(MAXL,MAXM), the temperature. ! ! Input/output, real TE(MAXL,MAXM), the turbulent epsilon. ! ! Input/output, real TK(MAXL,MAXM), the turbulent K. ! ! Input/output, real U(MAXL,MAXM), the horizontal velocity. ! ! Input/output, real V(MAXL,MAXM), the vertical velocity. ! ! Input/output, real W(MAXL,MAXM), the axial velocity. ! ! Input/output, real X(MAXL,MAXM), the X coordinates of the primary ! nodes. ! ! Input/output, real XC(MAXL,MAXM), the X coordinates of the corner ! nodes. ! ! Input/output, real Y(MAXL,MAXM), the Y coordinates of the primary ! nodes. ! ! Input/output, real YC(MAXL,MAXM), the Y coordinates of the corner ! nodes. ! implicit none ! integer maxbot integer maxl integer maxm ! real e(maxl,maxm) real gamt(maxl,maxm) integer i integer j integer jcrys1 integer jcrys2 integer l integer kcrys(maxl,maxm) integer kmelt(maxl,maxm) integer kvoid(maxl,maxm) integer m integer nbot logical ncflag(maxl,maxm) logical npflag(maxl,maxm) real p(maxl,maxm) real pc(maxl,maxm) real psi(maxl,maxm) real rueta(maxl,maxm) real ruksi(maxl,maxm) real t(maxl,maxm) real te(maxl,maxm) real tk(maxl,maxm) real u(maxl,maxm) real v(maxl,maxm) real vmag(maxl,maxm) real vort(maxl,maxm) logical vpflag(maxl,maxm) real w(maxl,maxm) real xbot(maxbot) real xc(maxl,maxm) real xp(maxl,maxm) real ybot(maxbot) real yc(maxl,maxm) real yp(maxl,maxm) ! m = (m+1)/2 call rhalf(e,l,m,maxl,maxm) call rhalf(gamt,l,m,maxl,maxm) call ihalf(kcrys,l,m,maxl,maxm) call ihalf(kmelt,l,m,maxl,maxm) call ihalf(kvoid,l,m,maxl,maxm) call lhalf(ncflag,l,m,maxl,maxm) call lhalf(npflag,l,m,maxl,maxm) call rhalf(p,l,m,maxl,maxm) call rhalf(pc,l,m,maxl,maxm) call rhalf(psi,l,m,maxl,maxm) call rhalf(rueta,l,m,maxl,maxm) call rhalf(ruksi,l,m,maxl,maxm) call rhalf(t,l,m,maxl,maxm) call rhalf(te,l,m,maxl,maxm) call rhalf(tk,l,m,maxl,maxm) call rhalf(v,l,m,maxl,maxm) call rhalf(u,l,m,maxl,maxm) call rhalf(vort,l,m,maxl,maxm) call lhalf(vpflag,l,m,maxl,maxm) call rhalf(w,l,m,maxl,maxm) call rhalf(xp,l,m,maxl,maxm) call rhalf(xc,l,m,maxl,maxm) call rhalf(yp,l,m,maxl,maxm) call rhalf(yc,l,m,maxl,maxm) ! ! Take care of data defining the extent of the crystal. ! jcrys1 = 1 jcrys2 = jcrys2-m+1 ! ! Take care of data defining the crucible shape. ! nbot = (nbot+1)/2 do i = 1,nbot xbot(i) = xbot(i+nbot-1) ybot(i) = ybot(i+nbot-1) end do ! ! Set the velocity magnitudes. ! do i = 1,l do j = 1,m vmag(i,j) = sqrt(u(i,j)**2+v(i,j)**2) end do end do return end subroutine hello ! !*****************************************************************************80 ! !! HELLO prints out the program name, date, and purpose. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! write ( *, * ) ' ' write ( *, * ) 'HELLO!' write ( *, * ) ' ' write ( *, * ) ' CPLOT' write ( *, * ) ' Last revised on 12 June 1996' write ( *, * ) ' ' write ( *, * ) ' CPLOT produces X Window, PostScript or CGM' write ( *, * ) ' graphics output from the data of the 2D finite ' write ( *, * ) ' volume program MASTRAPP2D.' write ( *, * ) ' ' write ( *, * ) ' Scalar quantities can be displayed using contour' write ( *, * ) ' lines or contour colors.' write ( *, * ) ' ' write ( *, * ) ' Vector quantities can be displayed as vectors,' write ( *, * ) ' vector direction fields, stream lines or stream' write ( *, * ) ' contours, and magnitude contours.' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) ' CHECK OUT THE NEW CRYSC command!' write ( *, * ) ' ' return end subroutine help ! !*****************************************************************************80 ! !! HELP prints out a list of commands. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! write ( *, * ) ' ' write ( *, * ) 'Commands:' write ( *, * ) ' ' write ( *, * ) 'AXis Show the axis of symmetry.' write ( *, * ) 'BO Show the boundary.' write ( *, * ) 'BAr Show color bar.' write ( *, * ) 'C Choose colors.' write ( *, * ) 'CC = Choose color table.' write ( *, * ) 'CP Corrected pressure contours.' write ( *, * ) 'CPC Corrected pressure colors.' write ( *, * ) 'CRUcible Show the crucible wall.' write ( *, * ) 'CRYB Show the crystal boundary.' write ( *, * ) 'CRYS Display items in crystal.' write ( *, * ) 'CRYSC Color in the crystal.' write ( *, * ) 'CV Show control volumes.' write ( *, * ) 'CVN Show control volume numbers.' write ( *, * ) 'DAT = Specify the plot data file.' write ( *, * ) 'DEV = Specify the graphics output device.' write ( *, * ) 'DIFF Show the difference between two files.' write ( *, * ) 'DOUBLE Reflect the region about the Y axis.' write ( *, * ) 'DXCDP Plot (dXCdP,dYCdP) for a parameter.' write ( *, * ) 'ECHO Echo user input to output file.' write ( *, * ) 'FILE = Set the output file name.' write ( *, * ) 'FRAME Show the frame.' write ( *, * ) 'FS Show the free surface.' write ( *, * ) 'G Show the current graph.' write ( *, * ) 'GRACE = Set the grace margin.' write ( *, * ) 'H Help (print this list).' write ( *, * ) 'HALF Don''t reflect the region about the Y axis.' write ( *, * ) 'ICMAX Set maximum available color.' write ( *, * ) 'ICMIN Set minimum available color.' write ( *, * ) 'JCMAX Set maximum used color.' write ( *, * ) 'JCMIN Set minimum used color.' write ( *, * ) 'MELT Display items in the melt.' write ( *, * ) 'MELTC Color in the melt (NOT IMPLEMENTED YET).' write ( *, * ) 'NB Show bottom nodes.' write ( *, * ) 'NC Show corner nodes.' write ( *, * ) 'NP Show primary nodes.' write ( *, * ) 'NCON = Set number of contour lines.' write ( *, * ) 'NOFRAME Do not show frame.' write ( *, * ) 'NXSKIP Set X node skip.' write ( *, * ) 'NYSKIP Set Y node skip.' write ( *, * ) 'OV Overlay next graphs.' write ( *, * ) 'P Pressure contours.' write ( *, * ) 'PC Pressure colors.' write ( *, * ) 'Q Quit.' write ( *, * ) 'S Stream lines.' write ( *, * ) 'SC Stream colors.' write ( *, * ) 'SCALECV = Set scale factor for control volume numbers.' write ( *, * ) 'SCALENB = Set scale factor for bottom nodes.' write ( *, * ) 'SCALENC = Set scale factor for corner nodes.' write ( *, * ) 'SCALENP = Set scale factor for primary nodes.' write ( *, * ) 'SCALEV = Set velocity vector scale.' write ( *, * ) 'TE Turbulent epsilon.' write ( *, * ) 'TEMP Temperatures.' write ( *, * ) 'TEMPC Temperature colors.' write ( *, * ) 'TITLE = Set title.' write ( *, * ) 'TITLE2 = Set subtitle.' write ( *, * ) 'TK Turbulent K.' write ( *, * ) 'U Unit velocities.' write ( *, * ) 'V Velocity vectors.' write ( *, * ) 'VIS Viscosity.' write ( *, * ) 'VM Velocity magnitude contours.' write ( *, * ) 'VMC Velocity magnitude colors.' write ( *, * ) 'VOID Display items in the void.' write ( *, * ) 'VOIDC Color in the void (NOT IMPLEMENTED YET).' write ( *, * ) 'VORT Vorticity contours.' write ( *, * ) 'VORTC Vorticity colors.' write ( *, * ) 'VPN Set visibile primary nodes.' write ( *, * ) 'W Angular momentum contours.' write ( *, * ) 'WC Angular momentum colors.' write ( *, * ) 'XC X coordinate contours.' write ( *, * ) 'XCC X coordinate colors.' write ( *, * ) 'YC Y coordinate contours.' write ( *, * ) 'YCC Y coordinate colors.' write ( *, * ) ' ' write ( *, * ) 'TH MH BH top, middle, bottom halves.' write ( *, * ) 'LH CH RH left, center, right halves.' write ( *, * ) ' ' write ( *, * ) 'TL TC TR top left, center, right quarters.' write ( *, * ) 'ML MC MR middle left, center, right quarters.' write ( *, * ) 'BL BC BR bottom left, center, right quarters.' write ( *, * ) 'X choose arbitrary X, Y subwindow.' write ( *, * ) ' ' write ( *, * ) 'FULL return to full picture.' write ( *, * ) ' ' return end subroutine iduble(ia,l,m,maxl,maxm) ! !*****************************************************************************80 ! !! IDUBLE shifts and copies integer data when the region is "doubled". ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer IA(MAXL,MAXM), the array to ! be adjusted. ! ! Input, integer L, M, the number of rows and columns of data. ! ! Input, integer MAXL, MAXM, the maximum allowed number of ! rows and columns of data. ! implicit none ! integer maxl integer maxm ! integer ia(maxl,maxm) integer i integer j integer jj integer l integer m ! ! Shift data to the right M-1 columns. ! do i = 1,l do jj = 2*m-1,m,-1 ia(i,jj) = ia(i,jj+1-m) end do end do ! ! Now "reflect" data into columns 1 to M-1. ! do i = 1,l do j = 1,m-1 ia(i,j) = ia(i,2*m-j) end do end do return end subroutine ihalf(ia,l,m,maxl,maxm) ! !*****************************************************************************80 ! !! IHALF adjusts an integer array when the region is cut in half. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer IA(MAXL,MAXM), the array to be shifted. ! ! Input, integer L, M, the number of rows and columns of data. ! ! Input, integer MAXL, MAXM, the maximum allowed number of ! rows and columns of data. ! implicit none ! integer maxl integer maxm ! integer ia(maxl,maxm) integer i integer j integer l integer m ! ! Shift data to the left M-1 columns. ! do i = 1,l do j = 1,m ia(i,j) = ia(i,j+m-1) end do end do return end subroutine init(b1jbl,b2jbl,b3jbl,cost,cost2,delx,dely,dev,echo,fildat, & filgrf,filinp,grace,icmax,icmin,icolor,icrys1,icrys2,iplot,itable,jcmax, & jcmin,jcrys1,jcrys2,l,lbar,m,maxbot,maxl,maxm,maxobj,nbot,ncflag,ncon, & nflag,npflag,nxskip,nyskip,object,ovrlay,p,pc,psi,reflec,scalecv,scalenb, & scalenc,scalenp,scalev,show,title,title2,u,v,vpflag,vmag,vort,x1max,x1min, & x2max,x2min,xbot,xc,xmax,xmin,xp,xsmax,xsmin,xtmax,xtmin,y1max,y1min,y2max, & y2min,ybot,yc,ymax,ymin,yp,ysmax,ysmin,ytmax,ytmin) ! !*****************************************************************************80 ! !! INIT initializes the values of data. ! ! ! Modified: ! ! 16 February 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! DEV Output, character*10 DEV. ! The graphics output device to be used. Current legal ! values for DEV include: ! ! cgmb - CGM binary file. ! ps - PostScript file. ! xws - X window screen (interactive). ! ! FILDAT Output, character*80 FILDAT. ! The name of the data file to be read in, which contains ! the information defining the mesh and the physical ! parameters. ! ! FILGRF Output, character*80 FILGRF, the name of the output ! graphics file. ! ! FILINP Output, character*80 FILINP. ! The name of a file in which will be placed a copy of the ! input typed by the user while running Display. ! ! GRACE Output, real GRACE. ! The size of the "grace" margin on the plot. ! ! ICMAX, ! ICMIN Output, integer ICMAX, ICMIN, the maximum and ! minimum color indices to use in the color bar. ! ! ICOLOR Output, integer ICOLOR(MAXOBJ). ! Contains the color indexes for each object. ! However, in some cases, ICOLOR is actual a color table ! index. ! ! ICRYS1 Output, integer ICRYS1, the I coordinate of the row of corner ! nodes that define the bottom of the crystal. ! ! ICRYS2 Output, integer ICRYS2, the I coordinate of the row of corner ! nodes that define the top of the crystal. ! ! IPLOT Output, integer IPLOT. ! The number of plots made so far. ! ! ITABLE Output, integer ITABLE, the desired color table. ! ! 1: low black to high white ! 2: low blue to high yellow ! 3: low red, high blue, with bands between. ! 4: low red, yellow, green, blue, high white. ! 5: low white, blue, green, yellow, high red. ! 6: low blue to high red ! 7: linear table between 2 user colors. ! 8: linear table between N user colors. ! 9: low white to high black. ! ! JCMAX, ! JCMIN Output, integer JCMAX, JCMIN, the maximum and ! minimum color indices to use in the color bar. ! ! LBAR Output, logical LBAR, is .TRUE. if the color bar should ! be shown. ! ! MAXBOT Input, integer MAXBOT, the maximum storage for the nodes ! that define the crucible bottom. ! ! MAXL, ! MAXM Input, integer MAXL, MAXM, the maximum allowed number of ! rows and columns of data. ! ! MAXOBJ Input, integer MAXOBJ. ! The number of graphical "objects". ! ! NBOT Output, integer NBOT, the number of crucible nodes. ! ! NCON Output, integer NCON. ! The number of contour lines to be drawn. This is ! initialized to 12, but may be changed by the user. ! ! NFLAG Output, logical NFLAG(MAXL,MAXM). ! ! NFLAG is used to "flag" which nodes are active, ! that is, to be displayed, and which not, in the graph. ! ! NXSKIP Output, integer NXSKIP. ! NXSKIP is used to "thin" out a vector plot. ! ! If NXSKIP = 1, then a standard vector plot is made. ! ! Otherwise, in the X direction, vectors are drawn only ! in columns 1, 1+NXSKIP, 1+2*NXSKIP and so on. ! ! NYSKIP Output, integer NYSKIP. ! NYSKIP is used to "thin" out a vector plot. ! ! If NYSKIP = 1, then a standard vector plot is made. ! ! Otherwise, in the Y direction, vectors are drawn only ! in rows 1, 1+NYSKIP, 1+2*NYSKIP and so on. ! ! OBJECT Input, character*25 OBJECT(MAXOBJ), the names of the ! graphical objects. ! ! OVRLAY Input, logical OVRLAY. ! If OVRLAY is true, then the next time that a plot is ! requested, a "new frame" command is suppressed, so that ! the new plot is shown on top of the previous one. ! ! P Output, real P(MAXL,MAXM), the pressure. ! ! PC Output, real PC(MAXL,MAXM), the corrected pressure. ! ! PSI Output, real PSI(MAXL,MAXM), the stream function. ! ! REFLEC Output, logical REFLEC, .TRUE. if the region is ! reflected, .FALSE. otherwise. ! ! SCALECV Output, real SCALECV. ! A scale factor for the control volume numbers, which has a ! default value of 1. ! ! SCALENB Output, real SCALENB. ! A scale factor for the size of the nodes that define the ! crucible shape, which has a default value of 1. ! ! SCALENC Output, real SCALENC. ! A scale factor for the size of the corner nodes, which has a ! default value of 1. ! ! SCALENP Output, real SCALENP. ! A scale factor for the size of the primary nodes, which has a ! default value of 1. ! ! SCALEV Output, real SCALEV. ! A scale factor for velocity vectors. This starts out at 1.0. ! ! SHOW Output, logical SHOW(MAXOBJ). ! Contains, for each object, a flag determining whether it ! is to be shown or not. ! ! TITLE Output, character*40 TITLE. ! A title for the plots. ! ! TITLE2 Output, character*40 TITLE2. ! A subtitle used in the profile plots. ! ! U Output, real U(MAXL,MAXM), the horizontal velocity. ! ! V Output, real V(MAXL,MAXM), the vertical velocity. ! ! VPFLAG Output, logical VPFLAG(MAXL,MAXM), a flag which determines ! whether vector data at a given node should be displayed. ! ! VMAG Output, real VMAG(MAXL,MAXM), the velocity magnitudes. ! ! X Output, real X(MAXL,MAXM), the X coordinates of the primary ! nodes. ! ! X1MAX, ! X1MIN Output, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! X2MAX, ! X2MIN Input, real X2MAX, X2MIN, the maximum and minimum X ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! XBOT Output, real XBOT(NBOT), the X coordinates of the crucible ! bottom nodes. ! ! XSMAX Input, real XSMAX. ! The maximum X coordinate of the data to be displayed. ! XSMAX defaults to XMAX, but can be made smaller to ! focus on a portion of the region. ! ! XSMIN Input, real XSMIN. ! The minimum X coordinate of the data to be displayed. ! XSMIN defaults to XMIN, but can be made larger to ! focus on a portion of the region. ! ! Y Output, real Y(MAXL,MAXM), the Y coordinates of the primary ! nodes. ! ! Y1MAX, ! Y1MIN Output, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! ! Y2MAX, ! Y2MIN Input, real Y2MAX, Y2MIN, the maximum and minimum Y ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! XBOT Output, real XBOT(NBOT), the X coordinates of the crucible ! bottom nodes. ! ! YSMAX Input, real YSMAX. ! The maximum Y coordinate of the data to be displayed. ! YSMAX defaults to YMAX, but can be made smaller to ! focus on a portion of the region. ! ! YSMIN Input, real YSMIN. ! The minimum Y coordinate of the data to be displayed. ! YSMIN defaults to YMIN, but can be made larger to ! focus on a portion of the region. ! implicit none ! integer maxbot integer maxl integer maxm integer maxobj ! real b1jbl(maxm) real b2jbl(maxm) real b3jbl(maxl) real cost real cost2 real delx real dely character ( len = 10 ) dev logical echo character ( len = 80 ) fildat character ( len = 80 ) filgrf character ( len = 80 ) filinp real grace integer i integer icmax integer icmin integer icolor(maxobj) integer icrys1 integer icrys2 integer iplot integer itable integer j integer jcmax integer jcmin integer jcrys1 integer jcrys2 integer l logical lbar integer m integer nbot logical ncflag(maxl,maxm) integer ncon logical nflag(maxl,maxm) logical npflag(maxl,maxm) integer nxskip integer nyskip character ( len = 25 ) object(maxobj) logical ovrlay real p(maxl,maxm) real pc(maxl,maxm) real psi(maxl,maxm) logical reflec real scalecv real scalenb real scalenc real scalenp real scalev logical show(maxobj) character ( len = 40 ) title character ( len = 40 ) title2 real u(maxl,maxm) real v(maxl,maxm) logical vpflag(maxl,maxm) real vmag(maxl,maxm) real vort(maxl,maxm) real x1max real x1min real x2max real x2min real xbot(maxbot) real xc(maxl,maxm) real xmax real xmin real xp(maxl,maxm) real xsmax real xsmin real xtmax real xtmin real y1max real y1min real y2max real y2min real ybot(maxbot) real yc(maxl,maxm) real ymax real ymin real yp(maxl,maxm) real ysmax real ysmin real ytmax real ytmi