! ! Graphics utility routines: ! module plt use input,only: icolor,igreyscale,ibox_clabs,iboxplt,ishadeneg, | mtime_del_mins implicit none ! ! input: real :: + pltmin=1.,pltmax=0., ! min,max of field to plot + conint=0., ! contour interval + scfac=1., ! scale factor + pltchsize=0. ! character size for axes labels integer :: + iblack=255, ! color index to black (255 -- see mkrgb) + iwhite=1, ! color index to white (1) + igks_cgm=0, ! gks workstation id for cgm + igks_ps=0, ! gks workstation id for ps + igks_x11=0, ! gks workstation id for X11 + iwk_cgm=0, ! ncarg workstation id for cgm + iwk_ps=0, ! ncarg workstation id for ps + iwk_x11=0, ! ncarg workstation id for X11 + iplt_cgm=0, ! current plot object id for cgm (HLU's) + iplt_ps=0 ! current plot object id for ps (HLU's) ! ! vectors on maps only: real :: + vlc=0., ! vector low cutoff + vhc=0., ! vector high cutoff + vsc=0., ! vector scale + vrl=0. ! vector realized length integer :: + labelv=0 ! label vector arrows ! ! output: real :: + zmin,zmax, ! min,max of field actually plotted + vmin,vmax, ! min,max of vector magnitudes + ciu ! contour interval used ! ! Work space for conpack, areas, etc: ! integer,parameter :: lrwrk=2000,liwrk=1000,liama=300000,lcra=10000 ! integer,parameter :: lrwrk=20000,liwrk=10000, ! | liama=300000,lcra=10000 ! integer,parameter :: lrwrk=200000,liwrk=100000, ! | liama=3000000,lcra=100000 integer,parameter :: lrwrk=2000000,liwrk=1000000, | liama=30000000,lcra=1000000 real :: rwrk(lrwrk),xcra(lcra),ycra(lcra) integer :: iwrk(liwrk),iama(liama),iarea(10),igrp(10) ! ! Color table: integer,parameter :: mxclrs=256 integer :: nclrs,ixclr(mxclrs) ! contains !------------------------------------------------------------------- subroutine setplt use proc,only: spval ! ! Set some plot parameters for the run: ! (this routine called only once, before any plots are made) ! call cpseti('LLP',3) ! line label positioning ! call cpseti('LLP',0) ! line label positioning (0 mean no line labels) call cpseti('LLO',1) ! line label orientation (1=local direction) call cpsetr('LLS',.015) ! line label size (default = .01) call cpsetc('HIT',' ') ! high label text string call cpsetc('LOT',' ') ! low label text string ! call cpseti('ILP',0) ! info label position flag (0=center of box) call cpsetc('ILT',' ') ! info label text string call mapstc('OU','CO') ! continental outlines call mapsti('DO',1) ! dotted outlines call cpsetr('ORV',1.e12) ! for polar and satv projection maps call vvsetr('USV -- U special value',spval) call vvsetr('VSV -- V special value',spval) call vvseti('SVF -- special value flag',3) call cpsetr('SPV -- special value',spval) call pcsetc('FC','$') return end subroutine setplt !------------------------------------------------------------------- subroutine contour(zz,idimx,nx,ny) use proc,only: spval ! ! Make 2-d contour: ! This routine does no labeling or axes outside plot window. ! Several parameters from pltopt.h are used ! ! Args: integer,intent(in) :: idimx,nx,ny real,intent(in) :: zz(idimx,ny) ! ! Locals: integer :: icls=12,i,j,ncl,iaia,iaib,ilev0,lty real :: clev,scalefac,zzz(idimx,ny),fmin,fmax logical :: is0lev ! true of there is a zero contour level real :: vl,vr,vb,vt,wl,wr,wb,wt ! write(6,"('Enter contour: idimx=',i4,' nx=',i4,' ny=',i4, ! | ' zz min,max=',2e12.4)") idimx,nx,ny,minval(zz), ! | maxval(zz) ! ! Put box around plot if requested: ! if (iboxplt > 0) call box(1) ! ! Add group id for label boxes: ! (Note label boxing does not work if shading negative areas, ! because label boxes come out solid black on hardcopy when ! using ctrans to convert cgm to ps) ! ! if (ibox_clabs > 0.or.(icolor <= 0.and.ishadeneg > 0)) ! + call cpseti('LLB',2) ! line label box (2=fill only) ! ! ! write(6,"('contour: ibox_clabs=',i2,' ishadeneg=',i2, ! + ' icolor=',i2)") ibox_clabs,ishadeneg,icolor ! ! if (ibox_clabs > 0) then ! if (icolor > 0.or.(icolor <= 0.and.ishadeneg > 0)) then ! call cpseti('GIL - group id for label boxes',5) ! call cpseti('LLB',2) ! line label box (2=fill only) ! endif ! endif if (ibox_clabs > 0) then if (icolor==0.and.ishadeneg>0) ibox_clabs = 0 endif if (ibox_clabs > 0) then call cpseti('GIL - group id for label boxes',5) call cpseti('LLB',2) ! line label box (2=fill only) endif ! ! Set scale factor if necessary: ! zzz = zz scalefac = scfac if (scalefac /= 1.) zzz = merge(zzz,zzz*scalefac,zzz==spval) ! ! Force cmin,cmax if necessary (also see iaib and iaia below if color) ! if (pltmin < pltmax) then call cpsetr('CMN -- contour minimum',pltmin) call cpsetr('CMX -- contour maximum',pltmax) else call cpsetr('CMN -- contour minimum',1.) call cpsetr('CMX -- contour maximum',0.) endif ! ! If conint <= 0., then let conpack pick contour levels (using icls as ! approx number of levels), otherwise force conint as contour level: ! (CLS (when > 0) is not used when CIS > 0., otherwise (CIS <= 0.) ! conpack will choose the interval such that at least CLS levels ! are used (here CLS == icls == 12)) ! call cpsetr('CIS -- contour interval specifier',conint) if (conint <= 0.) then call cpseti('CLS -- contour level selector',icls) else call cpseti('CLS -- contour level selector',1) endif ! ! Initialize conpack and pick contour levels: ! (ncl = number of contour levels) ! call getset(vl,vr,vb,vt,wl,wr,wb,wt,lty) call cprect(zzz,idimx,nx,ny,rwrk,lrwrk,iwrk,liwrk) call cppkcl(zzz,rwrk,iwrk) call cpgeti('NCL -- number of contour levels',ncl) ! ! Loop through contour levels: ! If contour level is < 0, make dashed contour lines unless we are ! shading negative areas (monochrome only), in which case draw ! solid black contour lines: ! ilev0 = 0 do i=1,ncl call cpseti('PAI -- parameter array index',i) call cpgetr('CLV -- contour level',clev) if (clev > -1.e-20 .and. clev < 1.e-20) ilev0 = i if (clev < -1.e-20) then if (ishadeneg > 0.and.icolor <= 0) then call cpseti('LLC -- line label color',iblack) call cpseti('CLC -- contour line color',iblack) call cpsetc('CLD -- contour line dash pattern', + '$$$$$$$$$$$$$$$$') else call cpsetc('CLD -- contour line dash pattern', + '$$$''''$$$''''') endif endif call cpseti('CLU -- contour level use flag',3) ! ! If shading, override default area identifiers: if (ishadeneg > 0.and.icolor <= 0) then call cpseti('AIA - area identifier above',0) call cpseti('AIB - area identifier below',0) endif ! ! If doing color fill, then always draw black contour lines, ! let line labels default white unless boxes are drawn, in ! which cases labels will be white over black-filled boxes. ! Also limit color fill to between pltmin and pltmax if necessary. ! (note mkrgb sets color index 255 to black for both cgm and ps) ! if (icolor > 0) then call cpseti('CLC -- contour line color',iblack) if (ibox_clabs.eq.0) then ! not boxing line labels call cpseti('LLC -- line label color',iblack) elseif (ibox_clabs.gt.0) then ! boxing line labels call cpseti('LLC -- line label color',iwhite) else ! if ibox_clabs < 0 then do not draw contour labels call cpseti('CLU -- contour level use flag',1) endif call cpseti('LLC -- line label color',iblack) call cpgeti('AIA - area identifier above',iaia) call cpgeti('AIB - area identifier below',iaib) if (pltmin < pltmax) then if (iaia.gt.ncl) call cpseti('AIA',-1) if (iaib.eq.1) call cpseti('AIB',-1) endif endif enddo ! i=1,ncl ! call arseti('RC',2) ! ncarg4.0 only ! ! Set up area maps, draw contour lines and labels: ! if (icolor <= 0) then ! monochrome call arinam(iama,liama) ! initialize area map ! ! If shading negative areas, add contour line at 0 if necessary, ! and set area id below 0 to 1 (see routine shader). If there was ! NOT already a 0 contour (ilev0 == 0), don't draw the line, otherwise ! draw the line double thickness: ! if (ishadeneg > 0) then ! shading areas < 0 if (ilev0 > 0) then ! ilev0 is already level 0. call cpseti('PAI - parameter array index',ilev0) call cpseti('CLU -- contour level use flag',3) ! call cpsetr('CLL -- contour level line width',2.) ! write(6,"('ilev0=',i2,' set CLU=3')") ilev0 else ! make new level 0. call cpseti('NCL - number of contour levels',ncl+1) call cpseti('PAI - parameter array index',ncl+1) call cpsetr('CLV - contour level value',0.) call cpseti('CLU -- contour level use flag',0) ! write(6,"('Made 0 level at ncl+1=',i2)") ncl+1 endif call cpseti('AIA - area identifier above',2) call cpseti('AIB - area identifier below',1) call cpclam(zzz,rwrk,iwrk,iama) ! add contour lines to am call arscam(iama,xcra,ycra,lcra,iarea,igrp,10,shader) endif call cplbam(zzz,rwrk,iwrk,iama) ! add contour labels to am call cpcldm(zzz,rwrk,iwrk,iama,drawcl) ! draw contour lines call cplbdr(zzz,rwrk,iwrk) ! draw contour labels else ! color nclrs = ncl+2 call mkrgb(nclrs,ixclr) ! make color table call arinam(iama,liama) ! initialize area map call cpclam(zzz,rwrk,iwrk,iama) ! add contour lines to am call arscam(iama,xcra,ycra,lcra,iarea, ! scan and color am + igrp,10,colram) call arinam(iama,liama) ! clear area map call cplbam(zzz,rwrk,iwrk,iama) ! add label boxes to am call cpcldm(zzz,rwrk,iwrk,iama,drawcl) ! draw contour lines ! write(6,"('contour before cplbdr: zzz min,max=',2e12.4)") ! | minval(zzz),maxval(zzz) call cplbdr(zzz,rwrk,iwrk) ! draw contour labels endif ! ! Define zmin, zmax, and ciu (used in plot labeling): call cpgetr('ZMN -- field minimum',zmin) call cpgetr('ZMX -- field maximum',zmax) call cpgetr('CIU -- contour interval used',ciu) return end subroutine contour !------------------------------------------------------------------- subroutine drawcl(xcs,ycs,ncs,iai,iag,nai) ! integer,intent(in) :: ncs,nai real,intent(in) :: xcs(ncs),ycs(ncs) integer,intent(in) :: iai(nai),iag(nai) integer :: i,idr ! idr = 1 do i=1,nai if (iai(i).lt.0) idr=0 enddo if (idr.ne.0) call curved(xcs,ycs,ncs) return end subroutine drawcl !------------------------------------------------------------------- subroutine colram(xcra,ycra,ncra,iarea,igrp,ngrps) ! ! Fill areas: ! xcra(ncra),ycra(ncra) = points defining a polygonal area ! iarea(ngrps) = area identifiers ! igrp(ngrps) = group identifiers (groups 1-4 are default) ! ! Edge group 1=ocean, 2=continents, 3=contour levels, 4=vertical strips ! ! Args: integer,intent(in) :: ngrps,ncra real,intent(in) :: xcra(ncra),ycra(ncra) integer,intent(in) :: iarea(ngrps),igrp(ngrps) ! ! Locals: integer :: ifll,i ! ifll=1 do i=1,ngrps if (iarea(i).lt.0) ifll=0 enddo if (ifll.ne.0) then ! fill the area ifll=0 do i=1,ngrps if (igrp(i).eq.3) ifll=iarea(i) ! area id of a contour level enddo if (ifll.gt.0.and.ifll.le.nclrs) then ! within the color tabl call gsfaci(ixclr(ifll+1)) call gfa(ncra-1,xcra,ycra) elseif (ifll.gt.nclrs) then call gsfaci(nclrs+1) call gfa(ncra-1,xcra,ycra) endif endif return end subroutine colram !------------------------------------------------------------------- subroutine shader(xcs,ycs,ncs,iai,iag,nai) ! ! Args: integer,intent(in) :: ncs,nai real,intent(in) :: xcs(ncs),ycs(ncs) integer,intent(in) :: iai(nai),iag(nai) ! ! Locals: integer,parameter :: lenwrk=5000 real :: dst(lenwrk) integer :: ind(lenwrk),ish,i ! ! Shade area identifiers that were set to 1 (see AIA,AIB in contour) ish = 0 do i=1,nai if (iag(i).eq.3.and.iai(i).eq.1) ish = 1 ! write(6,"('shader: nai=',i3,' i=',i3,' iag(i)=',i3, ! + ' iai(i)=',i3,' ish=',i2)") nai,i,iag(i),iai(i),ish enddo if (ish.ne.0) then ! ! Dot pattern: ! call sfsetr('spacing of fill lines',.004) call sfseti('angle of fill lines',0) call sfseti('dot-fill flag',1) call sfwrld(xcs,ycs,ncs-1,dst,lenwrk,ind,lenwrk) ! ! 45 deg hatch pattern: ! ! call sfseti('ANGLE',45) ! call sfsetr('SPACING',.006) ! call sfwrld(xcs,ycs,ncs-1,dst,lenwrk,ind,lenwrk) ! call sfseti('ANGLE',135) ! call sfnorm(xcs,ycs,ncs-1,dst,lenwrk,ind,lenwrk) endif return end subroutine shader !------------------------------------------------------------------- subroutine box(itic) ! ! Arg: integer,intent(in) :: itic ! ! Locals: real :: ticlen=.015 real :: vl,vr,vb,vt,wl,wr,wb,wt,xpos,xpos1,ypos,ypos1 integer :: lty,i,ii,iii ! ! Draw box around perimeter of viewport; ! If itic=1 -> add tic marks to scale axes in tenths; ! ! Save set and get into fractional coords: call getset(vl,vr,vb,vt,wl,wr,wb,wt,lty) call set(0.,1.,0.,1.,0.,1.,0.,1.,1) ! ! Draw the box: call line(0.,0.,1.,0.) call line(1.,0.,1.,1.) call line(1.,1.,0.,1.) call line(0.,1.,0.,0.) ! ! Add tic marks for reference: if (itic.gt.0) then do 100 i=1,4 do 101 ii=1,10 if (i.eq.1) then ! left axis xpos = 0. xpos1 = xpos+0.5*ticlen ypos = 0.1*float(ii-1) do 102 iii=1,9 ypos = ypos + .01 call line(xpos,ypos,xpos1,ypos) 102 continue xpos1 = xpos+ticlen ypos = 0.1*float(ii) ypos1 = ypos elseif (i.eq.2) then ! right axis xpos = 1. xpos1 = xpos-0.5*ticlen ypos = 0.1*float(ii-1) do 103 iii=1,9 ypos = ypos + .01 call line(xpos,ypos,xpos1,ypos) 103 continue xpos1 = xpos-ticlen ypos = 0.1*float(ii) ypos1 = ypos elseif (i.eq.3) then ! bottom axis ypos = 0. ypos1 = ypos+0.5*ticlen xpos = 0.1*float(ii-1) do 104 iii=1,9 xpos = xpos + .01 call line(xpos,ypos,xpos,ypos1) 104 continue ypos1 = ypos+ticlen xpos = 0.1*float(ii) xpos1 = xpos elseif (i.eq.4) then ! top axis ypos = 1. ypos1 = ypos-0.5*ticlen xpos = 0.1*float(ii-1) do 105 iii=1,9 xpos = xpos + .01 call line(xpos,ypos,xpos,ypos1) 105 continue xpos = 0.1*float(ii) xpos1 = xpos ypos1 = ypos-ticlen endif call line(xpos,ypos,xpos1,ypos1) 101 continue 100 continue endif call set(vl,vr,vb,vt,wl,wr,wb,wt,lty) return end subroutine box !------------------------------------------------------------------- subroutine mkrgb(numclr,iclrs) ! ! This subroutine calculates the red-green-blue mixtures based on ! equations in "Color Spaces for Computer Graphics" by Joblove and ! Greenberg from "Communication of the Association for Computing ! Machinery--Proceedings of SIGGRAPH, 1978". ! On input -- ! numclr ------- integer number of distinct colors for color bar ! if = -14 use the special color bar for transparencies ! if = 14 use the special color bar for slides ! On output -- ! The color table will be set up going from blues to pink. ! ! Args: integer,intent(inout) :: numclr integer,intent(out) :: iclrs(numclr) ! real rgb(3,mxclrs),rgbp(3,mxclrs),rgbt(3,16) real,parameter :: pi = 3.141592654 real :: one6th=1./6., + two6th=2./6., + thr6th=3./6., + for6th=4./6., + fiv6th=5./6., + twn436=24./36. ! color offset real :: hdif,h,r,g,b integer :: i,j,index ! ! RGB for 14 color bar (for slides) ! data ((rgb(i,j),i=1,3),j=1,10) / 1 0.0, 0.0, 0.0, 2 1.0, 1.0, 1.0, 3 0.0, 0.0, 0.9, 4 0., 0.1667, 1.0, 5 0.0, 0.5, 1.0, 6 0.0, 0.75, 1.0, 7 0.0, 0.5, 0.0, 8 0.0 , 0.80, 0.0, 9 0.33, 1.0, 0.0, 1 1.0, 1.0, 0.0/ data ((rgb(i,j),i=1,3),j=11,16) / 1 1.0, 0.60, 0.00, 2 1.0, 0.33, 0.00, 3 1.0, 0.0, 0.00, 4 1.0, 0.0, 0.70, 5 1.0, 0.0, 1.00, 6 1.0, 0.3, 1.00/ ! RGB for color transparencies data ((rgbt(i,j),i=1,3),j=1,10) / 1 0.0, 0.0, 0.0, 2 1.0, 1.0, 1.0, 3 0.0, 0.0, 1.0, 4 0.334, 0.667, 1.0, 5 0.0, 1.0, 1.0, 6 0.667, 1.0, 1.0, 7 0.334, 1.0, 0.667, 8 0.0, 1.0, 0.0, 9 0.667, 1.0, 0.0, 1 1.0, 1.0, 0.0/ data ((rgbt(i,j),i=1,3),j=11,16) / 1 1.0, 0.667, 0.334, 2 1.0, 0.334, 0.00, 3 1.0, 0.0, 0.334, 4 0.667, 0.0, 0.667, 5 1.0, 0.0, 1.0, 6 1.0, 0.667, 1.0/ !----------------------------------------------------------------------- ! check that the number of colors specified does not exceed the maximum. !----------------------------------------------------------------------- if (numclr.gt.mxclrs-2) then write(6,"('>>> mkrgb warning: number of colors requested = ', + i4,' but max allowed = ',i4,' will return max only')") + numclr,mxclrs numclr = mxclrs end if !----------------------------------------------------------------------- ! set fill area to solid !----------------------------------------------------------------------- call gsfais(1) !----------------------------------------------------------------------- ! calculate the rgb prime values. !----------------------------------------------------------------------- hdif=(34./36.)/float(numclr) h=twn436+(numclr/12.)*hdif do i=1,numclr h=h-hdif if (h.lt.0.) h=1. if (h.le.one6th) then rgbp(1,i)=1. rgbp(2,i)=6.*h rgbp(3,i)=0. elseif (h.le.two6th) then rgbp(1,i)=2.-6.*h rgbp(2,i)=1. rgbp(3,i)=0. elseif (h.le.thr6th) then rgbp(1,i)=0. rgbp(2,i)=1. rgbp(3,i)=6.*h-2. elseif (h.le.for6th) then rgbp(1,i)=0. rgbp(2,i)=4.-6.*h rgbp(3,i)=1. elseif (h.le.fiv6th) then rgbp(1,i)=6.*h-4. rgbp(2,i)=0. rgbp(3,i)=1. else rgbp(1,i)=1. rgbp(2,i)=0. rgbp(3,i)=6.-6.*h end if end do !----------------------------------------------------------------------- ! calculate the true rgb values using a sinusoidal interpolation !----------------------------------------------------------------------- do i=1,numclr do j=1,3 rgb(j,i)=( 1. + cos ( (1.-rgbp(j,i)) * pi) ) / 2. end do end do !----------------------------------------------------------------------- ! set the GKS number of colors to numclr and set the GKS color bar !----------------------------------------------------------------------- if (igks_cgm.gt.0) then ! write(6,"('mkrgb call gscr: igks_cgm=',i3,' iblack=',i3)") ! + igks_cgm,iblack call gscr(igks_cgm,iblack,0.,0.,0.) endif if (igks_ps.gt.0) call gscr(igks_ps,iblack,0.,0.,0.) do index = 2,numclr r = rgb(1,index) g = rgb(2,index) b = rgb(3,index) if(igreyscale > 0) then r = 0.9 - 0.8*(float(index)-2.)/(float(numclr)-1.) g = r b = r endif iclrs(index) = + 2 + ifix(252.*(float(index)-2.)/(float(numclr)-1.)) if (igks_cgm.gt.0) then ! write(6,"('mkrgb call gscr: index=',i3,' numclr=',i3, ! + ' iclrs(index)=',i3)") index,numclr,iclrs(index) call gscr(igks_cgm,iclrs(index),r,g,b) endif if (igks_ps.gt.0) + call gscr(igks_ps,iclrs(index),r,g,b) end do return end subroutine mkrgb !------------------------------------------------------------------- subroutine labrect(xx,nx,yy,ny,xlab,ylab,chsize) ! ! Draw axes and axis numeric and char labels for any rectangular ! plot (not called for elliptical projections) -- uses autograph ! ! On input: ! xx(nx) array of nx x-coord values (e.g., -180 -> 180) ! yy(ny) array of ny y-coord values (e.g., -87.5 -> 87.5) ! ! Args: integer,intent(in) :: nx,ny real,intent(in) :: xx(nx),yy(ny),chsize character(len=*),intent(in) :: xlab,ylab ! ! AGORIP common block obtained from D. Kennison 12/95. ! If the 3rd through 8th of these is set to 4 (rather than default of ! 8 in AG), then char sizes can be smaller (for multiplt): ! real :: smrl integer :: isld,mwcl,mwcm,mwce,mdla,mwcd,mwdq,inif common /agorip/ smrl,isld,mwcl,mwcm,mwce,mdla,mwcd,mwdq,inif ! ! Locals: integer :: length real :: charsz ! call agseti("SET.",4) if (chsize.le.0.) then charsz = .025 else charsz = chsize endif mwcl=4 mwcm=4 mwce=4 mdla=4 mwcd=4 mwdq=4 c write(6,"('labrect: set values in /agorip/ to 4 charsz=',f9.4)") c + charsz c c x-axis: c call agsetc("LABEL/NAME.","B") if (nx.gt.0) then call agsetr("AXIS/BOTTOM/CONTROL.",1.) call agseti("LINE/NUMBER.",-100) call agsetf("LINE/CHARACTER.",charsz) length = len(xlab) call agseti("LINE/MAX.",length) call agsetc("LINE/TEXT.",xlab) call agsetf ('AXIS/BOTTOM/NUMERIC/WIDTH/MANTISSA.', charsz) else call agsetr("LABEL/DEF/SUPPRESSION.",1.) call agsetr("AXIS/BOTTOM/CONTROL.",-1.) call agsetc("LABEL/NAME.","T") call agsetr("AXIS/TOP/CONTROL.",-1.) endif c c y-axis: c call agsetc("LABEL/NAME.","L") if (ny.gt.0) then call agseti("LINE/NUMBER.",100) call agsetf("LINE/CHARACTER.",charsz) length = len(ylab) call agseti("LINE/MAX.",length) call agsetc("LINE/TEXT.",ylab) call agsetf ('AXIS/LEFT/NUMERIC/WIDTH/MANTISSA.', charsz) else call agsetr("LABEL/DEF/SUPPRESSION.",1.) call agsetr("AXIS/LEFT/CONTROL.",-1.) call agsetc("LABEL/NAME.","R") call agsetr("AXIS/TOP/CONTROL.",-1.) endif c c Disable top label: c call agsetc("LABEL/NAME.","T") call agsetr("LABEL/DEF/SUPPRESSION.",1.) c c Draw the background c ! write(6,"('labrect before agstup: nx=',i3,' xx=',/(8f9.2))") ! + nx,xx ! write(6,"('labrect before agstup: ny=',i3,' yy=',/(8f9.2))") ! + ny,yy call agstup(xx,1,0,nx,1, yy,1,0,ny,1) call agback c c Reenable labels: c call agsetc("LABEL/NAME.","L") call agsetr("LABEL/DEF/SUPPRESSION.",0.) call agsetc("LABEL/NAME.","R") call agsetr("LABEL/DEF/SUPPRESSION.",0.) call agsetc("LABEL/NAME.","B") call agsetr("LABEL/DEF/SUPPRESSION.",0.) call agsetc("LABEL/NAME.","T") call agsetr("LABEL/DEF/SUPPRESSION.",0.) c return end subroutine labrect !------------------------------------------------------------------- subroutine labaxes(nx,numx,mnrx,labx,ny,numy,mnry,laby,charsize) ! ! Label axes with integers, converted to chars: ! ! Args: integer,intent(in) :: nx,ny integer,intent(in) :: numx(nx),numy(ny) character(len=*),intent(in) :: labx,laby integer,intent(in) :: mnrx,mnry real,intent(in) :: charsize ! ! Locals: character(len=8) chnum real :: chsize,vl,vr,vb,vt,wl,wr,wb,wt,vpw,yoffs,xoffs, + ticlen,distmjr,xpos,ypos,distmnr integer :: lty,ix,iiy,iix,iy ! ! chsize passed in already scaled by viewport width: ! chsize = .014 if (charsize.gt.0.) chsize = charsize call getset(vl,vr,vb,vt, wl,wr,wb,wt, lty) c c Scale offsets and tic length by viewport width: c vpw = vr-vl yoffs = .04*vpw xoffs = .03*vpw ticlen = .012*vpw c c Draw perimeter around user frame: c call line(wl,wb,wr,wb) call line(wr,wb,wr,wt) call line(wr,wt,wl,wt) call line(wl,wt,wl,wb) c c Get in fractional coord system: c call set (0.,1.,0.,1.,0.,1.,0.,1.,1) c c X-axis: if (nx.le.0) goto 500 xpos = vl ypos = vb - xoffs distmjr = (vr-vl)/float(nx-1) do ix=1,nx-1 call wrchnum(numx(ix),chnum) call plch(xpos,ypos,trim(chnum),chsize,0.,0.) call line(xpos,vb,xpos,vb+ticlen) call line(xpos,vt,xpos,vt-ticlen) do iix=1,mnrx distmnr = distmjr / float(mnrx) xpos = xpos + distmnr call line(xpos,vb,xpos,vb+0.5*ticlen) call line(xpos,vt,xpos,vt-0.5*ticlen) enddo enddo call wrchnum(numx(nx),chnum) call plch(xpos,ypos,trim(chnum),chsize,0.,0.) call plch(0.5*(vr+vl),vb-2.1*xoffs,trim(labx), + chsize,0.,0.) c c Y-axis: 500 if (ny.le.0) goto 550 xpos = vl - yoffs ypos = vb distmjr = (vt-vb)/float(ny-1) do iy=1,ny-1 call wrchnum(numy(iy),chnum) call plch(xpos,ypos,trim(chnum),chsize,0.,0.) call line(vl,ypos,vl+ticlen,ypos) call line(vr,ypos,vr-ticlen,ypos) do iiy=1,mnry distmnr = distmjr / float(mnry) ypos = ypos + distmnr call line(vl,ypos,vl+0.5*ticlen,ypos) call line(vr,ypos,vr-0.5*ticlen,ypos) enddo enddo call wrchnum(numy(ny),chnum) call plch(xpos,ypos,trim(chnum),chsize,0.,0.) call plch(vl-2.*yoffs,0.5*(vt+vb),trim(laby), + chsize,90.,0.) c 550 call set(vl,vr,vb,vt, wl,wr,wb,wt, lty) return end subroutine labaxes !------------------------------------------------------------------- subroutine wrchnum(num,ch) ! integer,intent(in) :: num character(len=*),intent(out) :: ch character(len=8) form integer :: min(8)=(/0,-9,-99,-999,-9999,-99999,-999999,-9999999/) integer :: max(8)=(/9,99,999,9999,99999,999999,9999999,99999999/) integer :: i ! do i=1,8 if (num.ge.min(i).and.num.le.max(i)) then write(form,"('(I',i1,')',4x)") i write(ch,form) num return endif enddo write(6,"('>>> wrchnum: num to big or too small =',i20)") num return end subroutine wrchnum !------------------------------------------------------------------- subroutine wrlab(lab,xpos,ypos,chsize,vp,ihq) ! ! Draw text label lab on plot centered at normalized coord xpos,ypos ! ! Args: character(len=*),intent(in) :: lab real,intent(in) :: xpos,ypos,chsize,vp(4) integer,intent(in) :: ihq ! ! Locals: real :: vl,vr,vb,vt,wl,wr,wb,wt,charsz integer :: lnlg ! call pcsetc('FC','$') call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) call set(0.,1.,0.,1.,0.,1.,0.,1.,1) if (chsize.le.0.) then charsz = .04 else charsz = chsize endif charsz = charsz*(vp(2)-vp(1)) if (ihq.gt.0) then call plch(xpos,ypos,lab,charsz,0.,0.) else call plch(xpos,ypos,lab,charsz,0.,0.) endif call set(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) return end subroutine wrlab !------------------------------------------------------------------- ! Sub wrlab6 is in wrlab6.f !------------------------------------------------------------------- subroutine sltxax(glonl,glonr,ut,yoffset) ! ! Add local time axis to existing longitude axis, at given ut ! Existing lon axis range is glonl (left) to glonr (right) ! (glonl,glonr given in the range -180. to 180.) ! This is usually called with glonl=glonr meaning 360 degrees longitude ! is being plotted, but it should work with narrower range of lons. ! ! Args: real,intent(in) :: glonl,glonr,ut,yoffset ! ! Locals: character(len=8) :: str=' ' character(len=32) :: xlab= 'LOCAL TIME (HRS)' integer :: i,ii,nmnr,idel,log,len integer :: nmjr=6 real :: yoff,chsize,ticlen,ynumoff,glon0,glon1,dlon,dhrs, + slt0,wlslt,slt,xpos0,xpos1,dmnr,yposax,ypos,vl,vr,vb,vt, + wl,wr,wb,wt,dum,xmnr,sltlab,xpos ! ! Externals: real,external :: fslt,rnd c c 12/95: c Yoffset from user is fraction of viewport height below bottom axis. c Chsize, ticlen, ynumoff are set as fractions of viewport width c call getset(vl,vr,vb,vt, wl,wr,wb,wt, log) yoff = yoffset*(vt-vb) chsize = .018*(vr-vl) ticlen = .02*(vr-vl) ynumoff = .05*(vr-vl) c glon0 = glonl+180. glon1 = glonr+180. dlon = glon1-glon0 if (dlon.eq.0.) dlon = 360. if (dlon.lt.0.) dlon = glon1+(360.-glon0) dhrs = dlon/15. ! delta hours across axis slt0 = fslt(dum,ut,glonl,1) wlslt = slt0 idel= int(dhrs)/nmjr ! delta hrs between major slt tics if (idel.le.0) idel = 1 nmnr = idel c c Below error fixed 12/13/95: c if (slt0-ifix(slt0).ne.0.) slt0 = rnd(slt0,del) if (slt0-ifix(slt0).ne.0.) slt0 = rnd(slt0,float(idel)) c call set(0.,1.,0.,1.,0.,1.,0.,1.,1) slt = slt0 xpos0 = vl + (slt-wlslt)*(vr-vl)/dhrs xpos1 = vl + (slt+idel-wlslt)*(vr-vl)/dhrs dmnr = 1./nmnr * (xpos1-xpos0) yposax = vb-yoff call line(vl,yposax,vr,yposax) ypos = vb-yoff-ynumoff do i=1,100 xpos = vl + (slt-wlslt)*(vr-vl)/dhrs if (xpos.lt.vl) goto 100 if (xpos.gt.vr) goto 101 call line(xpos,yposax,xpos,yposax+ticlen) sltlab = slt if (sltlab.ge.24.) sltlab = sltlab-24. write(str,"(f5.1)") sltlab call plch(xpos,ypos,trim(adjustl(str)),chsize,0.,0.) call line(xpos,yposax,xpos,yposax+ticlen) do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.ge.vl.and.xmnr.le.vr) + call line(xmnr,yposax,xmnr,yposax+.5*ticlen) enddo 100 continue slt = slt+idel enddo 101 continue call plch(0.5*(vl+vr),vb-yoff-2.0*ynumoff, + xlab(1:len_trim(xlab)),chsize,0.,0.) call set(vl,vr,vb,vt, wl,wr,wb,wt, log) return end subroutine sltxax !------------------------------------------------------------------- subroutine yaxzpht(hts1d,hts2d,gcmzp,id1,kmx,dzp,yy,ny,p0, + iht,iyax,vp,spval) ! ! Add extra non-linear y-axis on right in height or pressure: ! ! Args: integer,intent(in) :: id1,kmx,ny,iht,iyax real,intent(in) :: yy(ny),hts1d(kmx),hts2d(id1,kmx),gcmzp(kmx), + vp(4),p0,spval,dzp ! ! Locals: real :: yaxht(kmx),pmbht(ny),yaxmb(kmx) integer,parameter :: npresslab=14,nhtlab=6 real :: htlab(nhtlab),chsize,offtoax,offtolab,rndval,delrnd, + toprnd,botrnd,del integer :: izp,izp0,izp1,k,i,nzp,nminor,nnans ! ! Externals: integer,external :: ixfind real,external :: rnd ! ! For pressure (mb) right hand y-axis: note strings pfcode (with ! function codes) must match values of presslab (see yaxright.f) ! character(len=8) :: pfcode(npresslab) = + (/'500 ','300 ','10$S$2 ', + '10$S$1 ','10$S$0 ','10$S2$-1','10$S2$-2','10$S2$-3', + '10$S2$-4','10$S2$-5','10$S2$-6','10$S2$-7','10$S2$-8', + '10$S2$-9'/) real :: presslab(npresslab) = + (/ 500., 300., 1.e+2, + 1.e+1, 1.e+0, 1.e-1, 1.e-2, 1.e-3, + 1.e-4, 1.e-5, 1.e-6, 1.e-7, 1.e-8, + 1.e-9/) ! ! Check for nans: ! subroutine check_nans(f,id1,id2,id3,name,n_total,ispval,spval, ! iprint,ifatal) ! If ispval /= 0, then replace any nans w/ spval. ! Number of nans is returned in nnans. ! Return if any nans are found in the heights array: ! if (id1 > 0) then call check_nans(hts2d,id1,kmx,1,'ZSLICE',nnans, | 0,spval,0,0) else call check_nans(hts1d,kmx,1,1,'ZSLICE',nnans, | 0,spval,0,0) endif if (nnans > 0) then write(6,"('>>> WARNING yaxzpht: Found ',i8,' NaNs', | ' in hts array (out of ',i8,')')") | nnans,id1*kmx return endif ! call pcsetc('FC','$') ! ! If zp is on left y-axis (iht <= 0), then add non-linear ht y-axis ! yaxht(:) = 0. izp0 = 1 izp1 = kmx if (iht.le.0) then izp0 = ixfind(gcmzp,kmx,yy(1),dzp) izp1 = ixfind(gcmzp,kmx,yy(ny),dzp) endif do k=izp0,izp1 izp = k-izp0+1 if (id1.gt.0) then do i=1,id1 yaxht(izp) = yaxht(izp) + hts2d(i,k) enddo yaxht(izp) = yaxht(izp)/float(id1) else yaxht(izp) = hts1d(k) endif yaxmb(izp) = p0*exp(-gcmzp(k))*1.e-3 ! press(mb) enddo nzp = izp1-izp0+1 c c Offsets to new right y-axis from existing (unlabeled) right y-axis: c chsize = 0. offtoax = 0. offtolab = 0. c c Zp on lefthand y-axis, add righthand yaxes in height c (also add pressure(mb) axis if iyax > 1) c if (iht.le.0) then rndval = 10. if (yy(ny)-yy(1).le.5.) rndval = 5. botrnd = rnd(yaxht(1),rndval) toprnd = rnd(yaxht(ny),rndval) del = (toprnd-botrnd)/nhtlab delrnd = rnd(del,rndval) do k=1,nhtlab htlab(k) = botrnd + (k-1)*delrnd enddo nminor = ifix(delrnd)/10 if (delrnd.le.40.) nminor = ifix(delrnd)/5 if (delrnd.le.10.) nminor = 5 c c Change offsets if adding 2nd right hand y-axis: c call yaxright(yy,yaxht,ny,htlab,nhtlab,'(I4)',1,pfcode, + nminor,'AVE HEIGHT (KM)',offtoax,offtolab,chsize,vp,spval) c c Optionally add pressure(mb) axis in addition to height: c if (iyax.gt.1) then offtoax = .18 offtolab = offtoax+.14 nminor = 0 call yaxright(yy,yaxmb,ny,presslab,npresslab,' ',0,pfcode, + nminor,'PRESSURE (Mb)',offtoax,offtolab,chsize,vp,spval) endif c c Linear height is on lefthand y-axis, add righthand yaxis in pressure mb: c (use default offsets) c else call vecterp(yaxht,gcmzp,kmx,yy,pmbht,ny,spval,0) do k=1,ny if (pmbht(k).ne.spval) + pmbht(k) = p0*exp(-pmbht(k))*1.e-3 ! pressure(mb) at each height enddo nminor = 0 call yaxright(yy,pmbht,ny,presslab,npresslab,' ',0,pfcode, + nminor,'PRESSURE (Mb)',offtoax,offtolab,chsize,vp,spval) endif return end subroutine yaxzpht !------------------------------------------------------------------- subroutine yaxright(yax_l,yax_r,ny,ynum,nynum,form,int,fcode, + nmnr,lab,rofftoax,rofftolab,rchsize,vp,spval) c c Draw extra right-hand axis: c ! ! Args: integer,intent(in) :: ny,nmnr,nynum,int real,intent(in) :: yax_l(ny),yax_r(ny),ynum(nynum),vp(4),spval, + rofftoax,rofftolab,rchsize character(len=*),intent(in) :: form,lab,fcode(nynum) ! ! Locals: real :: vpw,ticlen,offtoax,offtonum,offtolab,chsize,vl,vr,vb,vt, + xpos,ypos,ynumcur,ynumprev,sl,sr,sb,st,wl,wr,wb,wt,del,sign, + yval integer :: k,nmjr,lty c c Offset from existing right axis to new right axis is given in normalized c coords. Set offsets to numeric labels and axis info label: c vpw = vp(2)-vp(1) ! viewport width ticlen = .012*vpw c offtoax = .02 if (rofftoax.gt.0.) offtoax = rofftoax offtoax = offtoax*vpw c offtonum = offtoax+.03*vpw offtonum = offtoax+.035*vpw c offtolab = .14 if (rofftolab.gt.0) offtolab = rofftolab offtolab = offtolab*vpw c chsize = .018 if (rchsize.gt.0.) chsize = rchsize chsize = chsize*vpw c c Get into normalized coords: c call getset(sl,sr,sb,st,wl,wr,wb,wt,lty) call set(0.,1.,0.,1.,0.,1.,0.,1.,1) vl = vp(1) vr = vp(2) vb = vp(3) vt = vp(4) c c Draw the axis: c xpos = vr+offtoax call line(xpos,vb,xpos,vt) c c Draw numeric labels, tics, and info label: c nmjr = 0 do k=1,nynum call vecterp(yax_r,yax_l,ny,ynum(k),yval,1,spval,0) if (yval.eq.spval) goto 100 ypos = vb+(vt-vb)*(yval-yax_l(1))/(yax_l(ny)-yax_l(1)) if (ypos.le.vt.and.ypos.ge.vb) then ynumcur = ynum(k) nmjr = nmjr+1 c write(6,"('k=',i2,' nmjr=',i2,' ynum(k)=',1pe9.3,' ypos=', c + f7.4,' yax_l(1),(ny)=',2(1pe9.2))") c + k,nmjr,ynum(k),ypos,yax_l(1),yax_l(ny) call line(xpos,ypos,xpos+ticlen,ypos) ! major tic mark call drnum(ynum(k),int,form,fcode(k), ! numeric label + vr+offtonum,ypos,chsize) c c Draw minor tic marks: c if (nmnr.gt.0.and.nmjr.ge.2) then c write(6,"('calling drmnr: ynum(k)=',1pe9.3,' ynumprev=', c + 1pe9.3)") ynum(k),ynumprev call drmnr(ynum(k),ynumprev,nmnr,vb,vt,yax_l,yax_r,ny, + xpos,ticlen,spval) del = ynum(k)-ynumprev if (nmjr.eq.2) then ! minor tics below bottom major sign = 1. if (del.gt.0.) sign = -1. call drmnr(ynumprev,ynumprev+sign*del,nmnr, + vb,vt,yax_l,yax_r,ny,xpos,ticlen,spval) endif endif ! nmjr >= 2 ynumprev = ynum(k) endif ! ypos for major is on axis 100 continue enddo c c Minor tics above top major: c if (nmnr.gt.0) then sign = 1. if (del.lt.0.) sign = -1. c write(6,"('above top major: ynumcur=',1pe9.2,' del=',1pe9.2, c + ' sign=',f6.2)") ynumcur,del,sign call drmnr(ynumcur+sign*del,ynumcur,nmnr, + vb,vt,yax_l,yax_r,ny,xpos,ticlen,spval) endif c c Axis string label: c call plch(vr+offtolab,0.5*(vt+vb),lab,chsize,90.,0.) c c Restore original set: c call set(sl,sr,sb,st,wl,wr,wb,wt,lty) return end subroutine yaxright !------------------------------------------------------------------- subroutine drnum(ynum,int,form,fcode,xpos,ypos,chsize) use input,only: label_quality ! character(len=*),intent(in) :: form,fcode integer,intent(in) :: int real,intent(in) :: xpos,ypos,chsize,ynum ! character*32 labnum character*8 label_quality_tmp if (len_trim(form).gt.0) then if (int.gt.0) then ! form is for integers write(labnum,form) ifix(ynum) else write(labnum,form) ynum endif call plch(xpos,ypos,trim(labnum),chsize,0.,-1.) else label_quality_tmp = label_quality label_quality = 'high' call plch(xpos,ypos,trim(fcode),chsize,0.,-1.) label_quality = label_quality_tmp endif return end subroutine drnum !------------------------------------------------------------------- subroutine drmnr(ytop,ybot,nmnr,vb,vt,yax_l,yax_r,ny,xpos, + ticlen,spval) ! ! Args: integer,intent(in) :: ny,nmnr real,intent(in) :: ytop,ybot,vb,vt,yax_l(ny),yax_r(ny), + xpos,ticlen,spval ! ! Locals: integer :: dmnr,i real :: ymnrval,ypos,yval c dmnr = (ytop-ybot)/float(nmnr) do i=1,nmnr-1 ! minor tics between pair of major tics if (dmnr.lt.0.) then ymnrval = ytop+float(i)*dmnr else ymnrval = ytop-float(i)*dmnr endif call vecterp(yax_r,yax_l,ny,ymnrval,yval,1,spval,0) if (yval.ne.spval) then ypos = vb+(vt-vb)*(yval-yax_l(1))/(yax_l(ny)-yax_l(1)) if (ypos.le.vt.and.ypos.ge.vb) then call line(xpos,ypos,xpos+0.5*ticlen,ypos) c write(6,"('drmnr on axis: i=',i2,' ytop,bot=', c + 2(1pe9.2),' dmnr=',1pe9.3,' ymnrval=',1pe9.2,' ypos=', c + f7.4)") i,ytop,ybot,dmnr,ymnrval,ypos else c write(6,"('drmnr off axis: i=',i2,' ytop,bot=', c + 2(1pe9.2),' dmnr=',1pe9.3,' ymnrval=',1pe9.2,' ypos=', c + f7.4)") i,ytop,ybot,dmnr,ymnrval,ypos endif else c write(6,"('dmnr yval==spval: i=',i2,' ymnrval=',1pe9.3, c + ' dmnr=',1pe9.3)") i,ymnrval,dmnr endif enddo return end subroutine drmnr !------------------------------------------------------------------- subroutine mkproj(proj,vp,cenlat,cenlon,rot,perimlat,eradii,icont, | pltlat0,pltlat1) ! ! Make an ezmap projection (type of projection in proj) ! ! On input: ! proj = 'CE' -> make cylindrical equidistant projection ! proj = 'ST' -> make stereographic polar projection ! proj = 'SV' -> make sat view projection ! proj = 'MO' -> make mollweide projection ! vp(4) = fractional coords left,right,bottom,top of ! desired location of projection ! perimlat = latitude of perimeter for stereographic proj ! pltlat0,pltlat1 = latitude range to plot on y-axis of CE projection ! ! Args: character(len=*),intent(in) :: proj real,intent(in) :: vp(4),cenlat,cenlon,rot,perimlat,eradii, | pltlat0,pltlat1 integer,intent(in) :: icont ! ! Locals: real :: vl,vr,vb,vt,wl,wr,wb,wt,ang real :: plm1(2),plm2(2),plm3(2),plm4(2) integer :: ity ! ! Cylindrical equidistant: ! if (proj.eq.'CE') then call mapint call mapsti('EL',0) call mapsti('PE',0) call mappos(vp(1),vp(2),vp(3),vp(4)) call mapset('MA',0.,0.,0.,0.) ! ! Namelist option fmap_latminmax not working at this time. ! It works if these lines are uncommented and the above mapset ! call ('MA') is not called, then it does work, but then ! map_censlt (i.e., cenlon) does NOT work. ! plm1=(/pltlat0,0./) ! plm2=(/-180,0./) ! plm3=(/pltlat1,0./) ! plm4=(/180.,0./) ! call mapset('CO',plm1,plm2,plm3,plm4) call maproj(proj,cenlat,cenlon,0.) ! note hardwired 0 rotation call mapint call getset(vl,vr,vb,vt,wl,wr,wb,wt,ity) call set(vp(1),vp(2),vp(3),vp(4),wl,wr,wb,wt,ity) ! ! Mollweide equidistant: ! elseif (proj.eq.'MO') then call mapint call mapsti('EL',1) call mapsti('PE',1) call mappos(vp(1),vp(2),vp(3),vp(4)) call mapset('MA',0.,0.,0.,0.) call maproj(proj,cenlat,cenlon,0.) call mapint call getset(vl,vr,vb,vt,wl,wr,wb,wt,ity) call set(vp(1),vp(2),vp(3),vp(4),wl,wr,wb,wt,ity) ! ! Stereographic (polar): ! elseif (proj.eq.'ST') then call mapsti('EL',1) call mappos(vp(1),vp(2),vp(3),vp(4)) call maproj(proj,cenlat,cenlon,rot) ang = 90.-abs(perimlat) call mapset('AN',ang,ang,ang,ang) call mapint ! ! Satellite view projection: ! elseif (proj.eq.'SV') then call mapsti('EL',1) call mapstc('OU','CO') call mapsti('DO',1) call mappos(vp(1),vp(2),vp(3),vp(4)) call mapstr('SA',eradii) call mapset('MA',0.,0.,0.,0.) call mapsti('LA',0) ! LA=0 suppresses labeling of meridians and poles call mapstr('GD',1.) ! distance (deg) between pts of grid lines ! call mapstr('GR',0.) ! GR=0 means no grid call mapstr('GR',20.) ! GR=20 means 20 deg grid spacing call mapint call maproj(proj,cenlat,cenlon,rot) call mapint ! ! Bad proj string: ! else write(6,"('>>> mkproj: unsupported projection: ',a)") proj endif if (icont.gt.0) call maplot return end subroutine mkproj !------------------------------------------------------------------- subroutine advframe(iwk1,igks1,iwk2,igks2,iwk3,igks3, + multiplt,iadvfr,nppf,msgout,type,iframe) ! ! Conditionally advance frame on each given workstation id that is > 0: ! (Typically, iwk1 is cgm, iwk2 is ps, iwk3 is x11) ! ! If multiplt <= 0, frame is advanced unconditionally. The message ! msgout is written to stdout. ! If multiplt > 0 and iadvfr <= 0 msgout and nppf are written to ! stdout (indicating the plot just put on the multiplt frame), ! but frame is *not* advanced. ! If multiplt > 0 and iadvfr > 0 msgout and nppf are written to stdout ! and frame is advanced. ! ! When frame is advanced, increment iframe. ! nppf = number of plots on current page. ! After a multiplt frame is advanced reset nppf to 0. ! type = plot type (written to stdout when advancing multiplt frame) ! ! When frame is to be advanced, only open workstations are advanced ! (iwk1 > 0 and/or iwk2 > 0). if the corresponding gks workstation ! is active the plotting was done with LLU's, and ngpict is used to ! advance the frame on the gks workstation. Otherwise the plotting ! was done with HLU's, and NhlFFrame is used to advance the ncarg ! workstation(s). ! ! Args: integer,intent(in) :: iwk1,igks1,iwk2,igks2,iwk3,igks3, + multiplt,iadvfr integer,intent(inout) :: iframe,nppf character(len=*),intent(in) :: msgout,type ! ! Locals: integer :: ier,iactive ! ! If not doing multiplt, write to stdout and advance frame: ! if (multiplt <= 0) then ! advance single-plot frame if (len_trim(msgout) > 0) then write(6,"('Frame ',i4,': ',a)") iframe,trim(msgout) else write(6,"('Frame ',i4,':')") iframe endif iframe = iframe+1 ! ! If doing multiplt, advance frame only if iadvfr > 0: ! else if (len_trim(msgout) > 0) + write(6,"('Frame ',i4,' (p',i2,'): ',a)") + iframe,nppf,trim(msgout) if (iadvfr > 0) then ! advance multi-plot frame write(6,"('Completed multi-plot frame ',i3,' (',i2, + ' plots) (',a,')')") iframe,nppf,trim(type) iframe = iframe+1 nppf = 0 else ! do not advance multiplt frame return endif endif ! ! Advance frame on one or more workstations: ! if (iwk1 > 0) then ! cgm metafile call gqwks(igks1,ier,iactive) if (iactive > 0) then ! for LLU's call ngpict(igks1,1) else ! for HLU's call NhlFFrame(iwk1,ier) endif endif if (iwk2 > 0) then ! ps file call gqwks(igks2,ier,iactive) if (iactive > 0) then ! for LLU's call ngpict(igks2,1) else ! for HLU's call NhlFFrame(iwk2,ier) endif endif if (iwk3 > 0) call NhlFFrame(iwk3,ier) ! X11 return end subroutine advframe !------------------------------------------------------------------- subroutine setmultivp(vp,iadvfr,iplot,multi,ipltrowcol,vpout) ! ! Given a viewport vp, return viewport vpout, adjusted according to ! iplot, multi, and ipltrowcol. Also return iadvfr=0 or 1 for whether ! or not frame is to be advanced (see sub advframe). ! ! Args: real,intent(in) :: vp(4) integer,intent(in) :: iplot,multi integer,intent(out) :: iadvfr real,intent(out) :: vpout(4) integer,intent(in) :: ipltrowcol(2) ! ! Locals: integer :: nrows,ncols,iframe,fncols,fnrows,irow,icol,i ! iadvfr = 0 if (multi.le.0) then do i=1,4 vpout(i) = vp(i) enddo return endif nrows = ipltrowcol(1) ncols = ipltrowcol(2) iframe = 0 fncols = float(ncols) fnrows = float(nrows) do irow = 1,nrows do icol = 1,ncols iframe = iframe+1 ! ! Calculate full frame viewport: ! vpout(1) = float(icol-1)/fncols vpout(2) = vpout(1)+1./fncols vpout(3) = float(nrows-irow)/fnrows vpout(4) = vpout(3)+1./fnrows ! ! Reduce to plot viewport: ! vpout(1) = vpout(1) + vp(1)/fncols vpout(2) = vpout(2) - (1.-vp(2))/fncols vpout(3) = vpout(3) + vp(3)/fnrows vpout(4) = vpout(4) - (1.-vp(4))/fnrows ! if (iframe.eq.iplot) then if (iframe.eq.nrows*ncols) iadvfr = 1 ! write(6,"('setmultivp: iplot=',i2,' iadvfr=', ! + i2,' vpout=',4f9.5)") iplot-1,iadvfr,vpout return endif enddo enddo write(6,"('>>> warning setmultivp: iplot=',i2)") iplot return end subroutine setmultivp !------------------------------------------------------------------- subroutine pltvec(u,v,id1,id2,x0,x1,y0,y1,incx,incy, + yboffset,iclr,iveclab,sclab,spv,map) ! ! Plot vectors on existing map projection: ! 2/97: add optional sclab to appear w/ scale label. ! ! Args: integer,intent(in) :: id1,id2,incx,incy,iveclab,iclr,map real,intent(in) :: u(id1,id2),v(id1,id2),x0,x1,y0,y1, + yboffset,spv character(len=*),intent(in) :: sclab ! ! Locals: integer :: ier,idm,iplci,ll,len real :: rdm,dmx,vl,vr,vb,vt,wl,wr,wb,wt ! call vvseti('MAP -- Mapping Flag',map) call vvseti('SET -- Set Call Flag',0) call vvsetr('XC1 -- Lower X Bound',x0) call vvsetr('XCM -- Upper X Bound',x1) call vvsetr('YC1 -- Lower Y Bound',y0) call vvsetr('YCN -- Upper Y Bound',y1) call vvsetr('USV -- U Special Value',spv) call vvsetr('VSV -- V Special Value',spv) call vvseti('SVF -- Special Value Flag',3) call vvseti('XIN -- X Axis Array Increment',incx) call vvseti('YIN -- Y Axis Array Increment',incy) call vvseti('VPO -- Vector Positioning Mode',1) call vvsetr('MXS -- Max vector text block char size',.015) call getset(vl,vr,vb,vt,wl,wr,wb,wt,ll) ! ! Label scale arrow and print vector stats if iveclab > 0: if (iveclab > 0) then call vvseti('LBL -- Vector label flag',1) call vvseti('VST -- Vector Statistics Output',1) endif ! ! Use vector package max text block and scale arrow (do not use ! min text block and scale arrow). ! Veclab is string to be printed with the scale arrow: ! (usually "UN+VN" or "UI+VI" to indicate neutrals or ions) ! (max length of MXT is 36) ! call vvsetc('MNT -- Min Vector Text String',' ') len = len_trim(sclab) if (len>0) call vvsetc('MXT -- Max Vector Text String', + trim(sclab(1:min(len,36)))) call vvseti('MXP -- Max Vector Text Block Positioning',3) call vvsetr('MXX -- Max Vector Text block x-coord',0.5) call vvsetr('MXY -- Max Vector Text block y-coord',-yboffset) ! ! Default arrow color index CTV is 0, meaning current gks polyline ! index will be used. Save current polyline index (iplci), and set ! it to iclr (which is white if plotting vectors only, or black if ! plotting vectors as overlay to color-fill contours). ! call gqplci(ier,iplci) call gsplci(iclr) ! ! Scale vector arrow and text always white: call vvseti('MNC -- min vector text block color',iwhite) call vvseti('MXC -- max vector text block color',iwhite) ! ! Set reference magnitude: ! (vsc is global to module plt. Is zero (i.e., is mag of vrl) unless ! vn_scale or vi_scale is set by user. ! call vvsetr('VRM -- Vector Reference Magnitude',vsc) ! ! Initialize vectors package, and get min,max magnitudes: ! (vmin,vmax are global to module plt and are used by mkmaplabs) ! call vvinit(u,id1,v,id1,rdm,idm,id1,id2,0.,0) call vvgetr('VMN -- Min Vector Magnitude',vmin) call vvgetr('VMX -- Max Vector Magnitude',vmax) ! ! Set magnitude and length of scale arrow: call vvgetr('DMX -- Get Device Max Vector Length',dmx) if (vrl.le.0.) vrl = 1.8*dmx/(vr-vl) call vvsetr('VRL -- Set Vector Realized Length',vrl) call vvsetr('VLC -- Set Vector Low Cutoff',-vlc) call vvsetr('VHC -- Set Vector High Cutoff',-vhc) ! ! Draw the vectors: call vvectr(u,v,rdm,idm,idm,rdm) call gsplci(iplci) ! restore gks polyline color index return end subroutine pltvec !------------------------------------------------------------------- subroutine labutxy(mtin,ntms,yy,ny,ylab,charsize,isltax,rlon) implicit none ! ! Label time x-axis; also label y-axis; ! Add local time axis below time axis if isltax > 0 ! Assumes constant delta time along x-axis ! ! Args: integer,intent(in) :: ntms,ny,mtin(3,ntms) integer,intent(inout) :: isltax character(len=*),intent(in) :: ylab real,intent(in) :: charsize,rlon,yy(ny) ! ! Locals: integer,parameter :: mxlabs=12 real,parameter :: rndmin = 5. real :: dday0,dday1,vl,vr,vb,vt,wl,wr,wb,wt,vpw,chsize, + dhrs,dmin,dum(1) integer :: log,ndays,nmnr,i,ii,n real :: dlabmin,dlabhrs,dlabday,dday integer :: mt(3,ntms) ! ! Externals: real,external :: cmt2dday,cmt2dhrs,cmt2dmin,rnd ! ! write(6,"('Enter labutxy: mtin=',/,(3i4))") mtin(:,1:ntms) mt = mtin ! init local mtimes do i=2,ntms if (mt(1,i) < mt(1,i-1)) then ! write(6,"('labutxy: i=',i3,' mt(1,i)=',i4,' mt(1,i-1)=', ! | i4,': changing mt(1,i) to ',i4)") i,mt(1,i), ! | mt(1,i-1),mt(1,i)+365 mt(1,i) = mt(1,i)+365 endif enddo dday0 = cmt2dday(mt(1,1),mt(2,1),mt(3,1)) dday1 = cmt2dday(mt(1,ntms),mt(2,ntms),mt(3,ntms)) ! write(6,"('labutxy: dday0=',f8.3,' dday1=',f8.3)") ! | dday0,dday1 c c Get set and scale chsize: c call getset(vl,vr,vb,vt, wl,wr,wb,wt, log) vpw = vr-vl chsize = .02 if (charsize.gt.0.) chsize = charsize chsize = chsize*vpw c c Let autograph label y-axis: c if (ny.gt.0) then call set(vl,vr,vb,vt,dday0,dday1,yy(1),yy(ny),1) call labrect(dum,0,yy,ny,' ',trim(adjustl(ylab)), + chsize/vpw+.01*vpw) else call set(vl,vr,vb,vt,dday0,dday1,0.,1.,1) endif ! ! Calculate delta time interval between ut labels: ! dhrs = cmt2dhrs(mt(1,ntms),mt(2,ntms),mt(3,ntms))- + cmt2dhrs(mt(1,1),mt(2,1),mt(3,1)) dmin = cmt2dmin(mt(1,ntms),mt(2,ntms),mt(3,ntms))- + cmt2dmin(mt(1,1),mt(2,1),mt(3,1)) dday = cmt2dday(mt(1,ntms),mt(2,ntms),mt(3,ntms))- + cmt2dday(mt(1,1),mt(2,1),mt(3,1)) ! write(6,"('labutxy: dhrs=',f12.3,' dmin=',f12.3,' dday=',f8.3)") ! | dhrs,dmin,dday dlabmin = dmin/float(mxlabs) dlabmin = rnd(dlabmin,rndmin) dlabhrs = dlabmin/60. dlabday = dlabmin/1440. nmnr = 1 ! ! subroutine labutax(mt,ntms,del,itime,nmnr,chsize) ! If itime = 1,2,3, -> del is in minutes, hours, days, or months ! ! Up to 4 hrs, label by minute (rounded to nearest minute): if (dhrs <= 4) then do i=3,10 if (mod(nint(dlabmin),i)==0) then nmnr = i exit endif enddo call labutax(mt,ntms,dlabmin,1,nmnr,chsize) ! ! 4 hrs to 4 days, label by hour (rounded to nearest hour): elseif (dhrs <= 4.*24.) then if (mod(ifix(dhrs),6)==0) dlabhrs = 6. dlabhrs = rnd(dlabhrs,1.) if (nint(dlabhrs)==1) then nmnr = 4 else do i=10,2,-1 if (mod(nint(dlabhrs),i)==0) then nmnr = i exit endif enddo endif call labutax(mt,ntms,dlabhrs,2,nmnr,chsize) ! ! 4 to 100 days, label by day (rounded to nearest day): elseif (dhrs <= 2400.) then dlabday = rnd(dlabday,1.) if (nint(dlabday)==1) then nmnr = 4 else do i=10,2,-1 if (mod(nint(dlabday),i)==0) then nmnr = i exit endif enddo endif call labutax(mt,ntms,dlabday,3,nmnr,chsize) ! ! 100 days to 6 months, label by day (rounded to nearest 5 days): elseif (dday <= 180.) then dlabday = rnd(dlabday,5.) nmnr = 5 do i=3,5 if (mod(nint(dlabday),i)==0) then nmnr = i exit endif enddo call labutax(mt,ntms,dlabday,3,nmnr,chsize) ! ! > 6 months, label by month (assume label every month, so ! delta month = 1). else nmnr = 4 call labutax(mt,ntms,1.,4,nmnr,chsize) endif ! ! slt axis: ! (draw slt axis only up to 10 days) ! if (isltax > 0) then if (dhrs.le.0.25) then call labsltax(rlon,.05,3,chsize) elseif (dhrs.le.0.75) then call labsltax(rlon,.25,5,chsize) elseif (dhrs.le.1.9) then call labsltax(rlon,.25,5,chsize) elseif (dhrs.le.5.0) then call labsltax(rlon,.75,3,chsize) elseif (dhrs.le.12.0) then call labsltax(rlon,2.,4,chsize) elseif (dhrs.le.24.0) then call labsltax(rlon,4.,4,chsize) elseif (dhrs.le.48.0) then call labsltax(rlon,6.,3,chsize) elseif (dhrs.le.72.0) then call labsltax(rlon,12.,6,chsize) elseif (dhrs.le.240.0) then call labsltax(rlon,24.,4,chsize) else isltax = 0 endif endif call set(vl,vr,vb,vt, wl,wr,wb,wt, log) end subroutine labutxy !------------------------------------------------------------------- subroutine labutax(mt,ntms,del,itime,nmnr,chsize) c c Plot numeric labels and tics on x-axis, at intervals of del time c If itime = 1,2,3,4 -> del is in minutes, hours, days, or months, c and labels should be in minutes, hours, days, or months. c (Also plot appropriate text label) c nmnr = number of minor spaces between each pair of labeled tics c charsize = character size (has default if 0 on input) c NOTE: on input, set call has been made, with wl,wr in decimal days ! ! Args: integer,intent(in) :: ntms,itime,nmnr,mt(3,ntms) real,intent(in) :: del,chsize ! ! Locals: real :: vl,vr,vb,vt,vpw,yoff,ticlen,ypos,rlab,xmnr,dmnr, | wl,wr,wb,wt,xpos,xpos0,xpos1,wlhrs,wrhrs,wlmin,wrmin, | xmon(12) integer :: monthday1(12) = ! first day of each month | (/1,32,60,91,121,152,182,213,244,274,305,335/) character(len=3) :: monthlab(12) = | (/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP', | 'OCT','NOV','DEC'/) integer :: idel,log,i,ii,md,mh,mm,mon,isfirst,imon,mday0,mday1, | monthday,iyear integer :: mt0(3),mt1(3) character(len=16) :: str character(len=120) :: mtimelab character(len=240) :: xlab ! ! Externals: real,external :: cmt2dhrs,cmt2dmin c call getset(vl,vr,vb,vt, wl,wr,wb,wt, log) call set(0.,1.,0.,1.,0.,1.,0.,1.,1) vpw = vr-vl c c Scale offsets to viewport (chsize was already scaled in labutxy): c yoff = .030*vpw ticlen = .02*vpw ypos = vb-yoff idel = int(del) c ! write(mtimelab,"(' (MTIMES ',i3.3,' ',i2.2,':',i2.2, ! + ' TO ',i3.3,' ',i2.2,':',i2.2,') ')") (mt(i,1),i=1,3), ! + (mt(i,ntms),i=1,3) ! ! 5/05: add delta time to label: mt0(:) = mt(:,1) mt1(:) = mt(:,ntms) if (mt1(1) > 366) mt1(1) = mt1(1)-366 write(mtimelab,"(' MTIMES ',i3.3,' ',i2.2,':',i2.2, | ' TO ',i3.3,' ',i2.2,':',i2.2,' BY ',i4,' MINS')") | mt0(:),mt1(:),mtime_del_mins c c Delta is in minutes: c if (itime.eq.1) then write(xlab,"('UT (MINUTES)',a)") trim(mtimelab) wlmin = cmt2dmin(mt(1,1),mt(2,1),mt(3,1)) wrmin = cmt2dmin(mt(1,ntms),mt(2,ntms),mt(3,ntms)) rlab = float(int(wlmin)/idel*idel) if (rlab.lt.wlmin) rlab = rlab+del xpos0 = vl + (rlab-wlmin)*(vr-vl)/(wrmin-wlmin) xpos1 = vl + (rlab+del-wlmin)*(vr-vl)/(wrmin-wlmin) dmnr = 1./nmnr * (xpos1-xpos0) do i=1,1000 xpos = vl + (rlab-wlmin)*(vr-vl)/(wrmin-wlmin) if (xpos.lt.vl.or.xpos.gt.vr) goto 300 do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.gt.vl) then call line(xmnr,vb,xmnr,vb+.5*ticlen) call line(xmnr,vt,xmnr,vt-.5*ticlen) endif enddo str = ' ' call dmintomt(rlab,md,mh,mm) write(str,"(i2)") mm str = trim(adjustl(str)) c write(6,"('i=',i2,' rlab=',f9.4,' str=',a,' xpos=',f9.4)") c + i,rlab,str,xpos call plch(xpos,ypos,trim(str),chsize,0.,0.) call line(xpos,vb,xpos,vb+ticlen) call line(xpos,vt,xpos,vt-ticlen) rlab = rlab+del c xposp = xpos enddo 300 continue do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.ge.vl.and.xmnr.le.vr) then call line(xmnr,vb,xmnr,vb+.5*ticlen) call line(xmnr,vt,xmnr,vt-.5*ticlen) endif enddo c c Delta is in hours: c elseif (itime.eq.2) then write(xlab,"('UT (HOURS)',a)") trim(mtimelab) wlhrs = cmt2dhrs(mt(1,1),mt(2,1),mt(3,1)) wrhrs = cmt2dhrs(mt(1,ntms),mt(2,ntms),mt(3,ntms)) rlab = float(int(wlhrs)/idel*idel) if (rlab.lt.wlhrs) rlab = rlab+del xpos0 = vl + (rlab-wlhrs)*(vr-vl)/(wrhrs-wlhrs) xpos1 = vl + (rlab+del-wlhrs)*(vr-vl)/(wrhrs-wlhrs) dmnr = 1./nmnr * (xpos1-xpos0) do i=1,1000 xpos = vl + (rlab-wlhrs)*(vr-vl)/(wrhrs-wlhrs) if (xpos.lt.vl.or.xpos.gt.vr) goto 200 do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.gt.vl) then call line(xmnr,vb,xmnr,vb+.5*ticlen) call line(xmnr,vt,xmnr,vt-.5*ticlen) endif enddo str = ' ' call dhrstomt(rlab,md,mh,mm) write(str,"(i2)") mh str = trim(adjustl(str)) c write(6,"('i=',i2,' rlab=',f9.4,' str=',a,' xpos=',f9.4)") c + i,rlab,str,xpos call plch(xpos,ypos,trim(str),chsize,0.,0.) call line(xpos,vb,xpos,vb+ticlen) call line(xpos,vt,xpos,vt-ticlen) rlab = rlab+del c xposp = xpos enddo 200 continue do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.ge.vl.and.xmnr.le.vr) then call line(xmnr,vb,xmnr,vb+.5*ticlen) call line(xmnr,vt,xmnr,vt-.5*ticlen) endif enddo c c Delta is in days c elseif (itime.eq.3) then xlab = ' ' write(xlab,"('UT (DAYS)',a)") trim(mtimelab) rlab = float(int(wl)/idel*idel) if (rlab.lt.wl) rlab = rlab+del xpos0 = vl + (rlab-wl)*(vr-vl)/(wr-wl) xpos1 = vl + (rlab+del-wl)*(vr-vl)/(wr-wl) dmnr = 1./nmnr * (xpos1-xpos0) do i=1,1000 xpos = vl + (rlab-wl)*(vr-vl)/(wr-wl) if (xpos.lt.vl.or.xpos.gt.vr) goto 100 do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.gt.vl) then call line(xmnr,vb,xmnr,vb+.5*ticlen) call line(xmnr,vt,xmnr,vt-.5*ticlen) endif enddo str = ' ' call ddaytomt(rlab,md,mh,mm) if (md.le.9) then write(str,"(i2)") md else write(str,"(i3)") md endif str = trim(adjustl(str)) call plch(xpos,ypos,trim(str),chsize,0.,0.) call line(xpos,vb,xpos,vb+ticlen) call line(xpos,vt,xpos,vt-ticlen) rlab = rlab+del c xposp = xpos enddo 100 continue do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.ge.vl.and.xmnr.le.vr) then call line(xmnr,vb,xmnr,vb+.5*ticlen) call line(xmnr,vt,xmnr,vt-.5*ticlen) endif enddo c c Delta is in months: ! Label beginning calendar day of each month applicable to the axis. c (assume model day == julian day) ! 10/26/09 btf: allow crossing year-boundary c elseif (itime.eq.4) then xlab = ' ' write(xlab,"('MONTH',a)") trim(mtimelab) xmon = 999. ! array op (12) mday0 = mt(1,1) mday1 = mt(1,ntms) imon = 0 iyear = 0 do ii=1,ntms if (mt(1,ii) == 366) iyear = iyear+1 do i=1,12 monthday = monthday1(i)+iyear*365 if (monthday == mt(1,ii)) then imon = imon+1 xmon(imon) = (vr-vl)*float(monthday-mday0) / | float(mday1-mday0) + vl if (xmon(imon) >= vl .and. xmon(imon) <= vr) then call plch(xmon(imon),ypos,monthlab(i),chsize,0.,0.) call line(xmon(imon),vb,xmon(imon),vb+ticlen) call line(xmon(imon),vt,xmon(imon),vt-ticlen) endif endif enddo enddo ! ! Add minor tics every .25 month: do i=1,11 if (xmon(i) /= 999..and.xmon(i+1) /= 999.) then do ii=1,3 xpos = xmon(i)+0.25*float(ii)*(xmon(i+1)-xmon(i)) call line(xpos,vb,xpos,vb+.5*ticlen) call line(xpos,vt,xpos,vt-.5*ticlen) if (i==1) then ! back to vl from first month xpos = xmon(i)-0.25*float(ii)*(xmon(i+1)-xmon(i)) if (xpos > vl) then call line(xpos,vb,xpos,vb+.5*ticlen) call line(xpos,vt,xpos,vt-.5*ticlen) endif endif if (i < 11) then ! forward to vr from last month if (xmon(i+2)==999.) then xpos = xmon(i+1)+0.25*float(ii)*(xmon(i+1)-xmon(i)) if (xpos < vr) then call line(xpos,vb,xpos,vb+.5*ticlen) call line(xpos,vt,xpos,vt-.5*ticlen) endif endif endif enddo endif enddo ! minor tics between months endif ! labelling by min,hr,day,month c c Add time label: c call plch(0.5*(vl+vr),vb-2.4*yoff,trim(xlab),chsize,0.,0.) call set(vl,vr,vb,vt, wl,wr,wb,wt, log) c return end subroutine labutax !------------------------------------------------------------------- subroutine labsltax(rlon,del,nmnr,chsize) c c Draw slt axis (hrs) below existing ut axis c (set call has been done with wl,wr = decimal days) c ! ! Args: real,intent(in) :: rlon,del,chsize integer,intent(in) :: nmnr ! ! Locals: real :: vl,vr,vb,vt,wl,wr,wb,wt,vpw,yoff,ticlen,ynumoff,ypos, + yposax,wlhrs,wlslt,wrhrs,wrslt,sltdhrs,slt,dmnr,xpos0,xpos1, + xpos,xmnr integer :: log,idel,i,ii,mdr,mhr,mmr,mml,mdl,mhl character(len=8) :: str character(len=32) :: xlab ! ! Externals: real,external :: cmt2dhrs,dhrs2ut,dut2dslt,rnd ! call getset(vl,vr,vb,vt, wl,wr,wb,wt, log) call set(0.,1.,0.,1.,0.,1.,0.,1.,1) vpw = vr-vl c c Scale offsets by viewport (chsize was already scaled by latutxy) c yoff = .14*vpw ticlen = .02*vpw ynumoff = .035*vpw c yposax = vb-yoff idel = int(del) if (idel.le.0) idel = 1 xlab = ' ' write(xlab,"('LOCAL TIME (HRS) (LON=',f7.2,')')") rlon c c Draw the axis: c call line(vl,yposax,vr,yposax) c c Convert to decimal slt (including model day): c call ddaytomt(wl,mdl,mhl,mml) wlhrs = cmt2dhrs(mdl,mhl,mml) wlslt = dut2dslt(wlhrs,rlon) call ddaytomt(wr,mdr,mhr,mmr) wrhrs = cmt2dhrs(mdr,mhr,mmr) wrslt = wlslt + (wrhrs-wlhrs) sltdhrs = float(int(wlslt)) sltdhrs = rnd(sltdhrs,del) slt = dhrs2ut(sltdhrs) if (sltdhrs.lt.wlslt) then do i=1,100 sltdhrs = sltdhrs+del if (sltdhrs.ge.wlslt) goto 50 enddo endif 50 continue ypos = vb-yoff-ynumoff xpos0 = vl + (sltdhrs-wlslt)*(vr-vl)/(wrslt-wlslt) xpos1 = vl + (sltdhrs+del-wlslt)*(vr-vl)/(wrslt-wlslt) dmnr = 1./nmnr * (xpos1-xpos0) c c Loop on labeled intervals: c do i=1,20 xpos = vl + (sltdhrs-wlslt)*(vr-vl)/(wrslt-wlslt) if (xpos.gt.vr) goto 100 if (xpos.ge.vl) then do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.ge.vl.and.xmnr.le.vr) + call line(xmnr,yposax,xmnr,yposax+.5*ticlen) enddo slt = dhrs2ut(sltdhrs) str = ' ' write(str,"(f5.2)") slt str = trim(adjustl(str)) c write(6,"('i=',i2,' slt=',f9.4,' sltdhrs=',f12.5,' str=',a, c + ' xpos=',f9.4)") i,slt,sltdhrs,str,xpos call plch(xpos,ypos,trim(str),chsize,0.,0.) call line(xpos,yposax,xpos,yposax+ticlen) endif sltdhrs = sltdhrs+del enddo 100 continue do ii=1,nmnr-1 xmnr = xpos-ii*dmnr if (xmnr.ge.vl.and.xmnr.le.vr) + call line(xmnr,yposax,xmnr,yposax+.5*ticlen) enddo call plch(0.5*(vl+vr),vb-yoff-2.2*ynumoff,xlab,chsize,0.,0.) call set(vl,vr,vb,vt, wl,wr,wb,wt, log) return end subroutine labsltax !------------------------------------------------------------------- subroutine pltvec1d(u,v,xx,nx,vscalemag,vscmag,yboffset,vp, + spval) ! ! Args: integer,intent(in) :: nx real,intent(in) :: u(nx),v(nx),xx(nx),vp(4),vscalemag,yboffset, + spval real,intent(out) :: vscmag ! ! Locals: real :: yy(2),vsum(nx,1),rdm,xmid,vrl, + vl,vr,vb,vt,wl,wr,wb,wt character(len=36) :: sclab = ' ' integer :: idm,ll,len ! ! Externals: real,external :: rnd interface function vecsum(u,v,id1,id2,spv) implicit none integer,intent(in) :: id1,id2 real,intent(in) :: u(id1,id2),v(id1,id2) real:: vecsum(id1,id2) ! result variable real,intent(in),optional :: spv end function vecsum end interface ! ! If scale magnitude not given, use maximum vector sum rounded ! to nearest 50. Y-axis will go from -vscmag to +vscmag. ! vscmag = vscalemag if (vscmag <= 0.) then vsum = vecsum(u,v,nx,1,spval) vscmag = rnd(maxval(vsum),50.) endif yy(1) = -vscmag yy(2) = vscmag idm = 0 rdm = 0. call set(vp(1),vp(2),vp(3),vp(4),xx(1),xx(nx),yy(1),yy(2),1) call vvseti('MAP -- Mapping Flag',0) call vvseti('SET -- Set Call Flag',0) call vvsetr('XC1 -- Lower X Bound',xx(1)) call vvsetr('XCM -- Upper X Bound',xx(nx)) call vvsetr('YC1 -- Lower Y Bound',0.) call vvsetr('YCN -- Upper Y Bound',0.) call vvseti('VPO -- Vector Positioning Mode',1) call vvseti('TRT -- Transformation Type',0) call vvsetc('MNT -- Min Vector Text String',' ') call vvsetc('MXT -- Max Vector Text String',' ') call vvseti('LBL -- Vector label flag',labelv) call vvsetr('LBS -- Vector label char size',.01) call vvsetr('AMX -- Arrow Head Max size',.02) call vvsetr('MXS -- Max vector text block char size',.015) ! ! Use vector package max text block and scale arrow (do not use ! min text block and scale arrow). ! Veclab is string to be printed with the scale arrow: ! (usually "UN+VN" or "UI+VI" to indicate neutrals or ions) ! (max length of MXT is 36) ! call vvsetc('MNT -- Min Vector Text String',' ') len = len_trim(sclab) if (len>0) + call vvsetc('MXT -- Max Vector Text String', + 'MAX '//trim(sclab(1:min(len,32)))) call vvseti('MXP -- Max Vector Text Block Positioning',3) call vvsetr('MXX -- Max Vector text block x-coord',0.5) call vvsetr('MXY -- Max Vector text block y-coord',-yboffset) call vvseti('MNC -- min vector text block color',iwhite) call vvseti('MXC -- max vector text block color',iwhite) call vvsetr('VRM -- Vector Reference Magnitude',vsc) call vvsetr('USV -- U special value',spval) call vvsetr('VSV -- V special value',spval) call vvseti('XIN -- X Axis Array Increment',1) call vvseti('YIN -- Y Axis Array Increment',1) ! call vvinit(u,nx,v,nx,rdm,idm,nx,1,rdm,idm) call vvgetr('VMN -- Min Vector Magnitude',vmin) call vvgetr('VMX -- Min Vector Magnitude',vmax) call getset(vl,vr,vb,vt,wl,wr,wb,wt,ll) xmid = 0.5*(vl+vr) ! ! Scale the arrows to match y-axis: ! (vrl is ndc length of max vector as frac of viewport width) vrl = vmax/vscmag * 0.5*(vt-vb)/(vr-vl) call vvsetr('VRL -- Set Vector Realized Length',vrl) call vvsetr('VLC -- Set Vector Low Cutoff',-vlc) call vvsetr('VHC -- Set Vector High Cutoff',-vhc) call vvectr(u,v,rdm,idm,idm,rdm) ! ! Draw horizontal line across middle of plot (arrow tails will ! lie on this line) call line(xx(1),0.,xx(nx),0.) end subroutine pltvec1d !------------------------------------------------------------------- subroutine pltconrec(plt,idx,nx,ny,xx,yy,xlab,ylab,iht, + flo,fhi,fint,vp,wl,wr,wb,wt,spv,iyax,fht,zp,nzp,p0) c c Contour using old conrec (tgcmproc calls this to plot phases, c using Cicely's modified conrec): c c /conre1/ is in conrec.a (old modified conrec) c /trans/ is in fx,fy below, (fx and fy are used only by velvct c (via pltvec) and conrec; conpack does not use fx and fy) c fx and fy (see below) are modified to allow iqzqz=2, so that they c are mapped to xmin,xmax,ymin,ymax and nset=1 is passed to conrec so c my set call is used. This way, I can do set call from -90 to +90 on c x-axis, but data will be plotted at -87.5 to +87.5. c ! ! Args: integer,intent(in) :: idx,nx,ny,iyax,nzp,iht real,intent(in) :: plt(idx,ny),xx(nx),yy(ny),vp(4), + fht(nx,nzp),zp(nzp),flo,fhi,fint,wl,wr,wb,wt,spv,p0 character(len=*),intent(in) :: xlab,ylab ! ! Commons (in conrec.a): integer :: ioffp,mm,nn,if,iqzqz real :: spval,xmin,xmax,ymin,ymax common /conre1/ ioffp,spval common/trans/mm,nn,xmin,xmax,ymin,ymax,if,iqzqz ! ! Locals: integer :: iyaxright real :: dzp,fmin,fmax c iqzqz = 2 mm = nx nn = ny xmin = xx(1) xmax = xx(nx) ymin = yy(1) ymax = yy(ny) ioffp = 1 spval = spv call set(vp(1),vp(2),vp(3),vp(4),wl,wr,wb,wt,1) call fminmax(plt,nx*ny,fmin,fmax,spval) zmin = fmin zmax = fmax ciu = fint ! IBM does not like -1430B arg to conrec. See comments re NDOT ! in conrec.f. ! call conrec(plt,idx,nx,ny,flo,fhi,fint,1,-1,-1430B,2) call conrec(plt,idx,nx,ny,flo,fhi,fint,1,-1,0,2) call labrect(xx,nx,yy,ny,xlab,ylab,pltchsize) iyaxright = iyax-iyax/10*10 ! ! Add extra right-hand axes in height and/or pmb: ! if (iyaxright.gt.0) then dzp = zp(2)-zp(1) call yaxzpht(zp,fht,zp,nx,nzp,dzp,yy,ny,p0,iht,iyaxright, + vp,spval) endif return end subroutine pltconrec !-------------------------------------------------------------------- subroutine pltstrm(f,fht,zp,xx,yy,nx,ny,nzp,dzp,glon,iyax,iht, + p0,vp,spval) ! ! Contour stream function: ! ! Args: integer,intent(in) :: nx,ny,nzp,iyax,iht real,intent(in) :: f(nx,ny),xx(nx),yy(ny),vp(4),fht(nx,nzp), + zp(nzp),p0,spval,glon,dzp ! ! Locals: integer,parameter :: nclev=13 real :: clev(nclev) = (/-1.E8,-1.E7,-1.E6,-1.E5,-1.E4,-1.E3, 0., + 1.E3, 1.E4, 1.E5, 1.E6, 1.E7, 1.E8/) real :: xmid,xleft,xright integer :: i,iyaxright ! ! Begin exec: xmid = 0.5*(vp(1)+vp(2)) call cpseti('SET',0) call cpseti('MAP',0) call cpsetr('XC1',xx(1)) call cpsetr('XCM',xx(nx)) call cpsetr('ILX',xmid) call cpsetr('ILY',-.16) call cpseti('ILP',0) call cpsetr('ILS',.016) call cpsetr('YC1',yy(1)) call cpsetr('YCN',yy(ny)) xleft = xx(1)-(xx(2)-xx(1))/2. xright = xx(nx)+(xx(2)-xx(1))/2. call set(vp(1),vp(2),vp(3),vp(4),xleft,xright,yy(1),yy(ny),1) call cpseti('CLS',0) call cpseti('NCL',nclev) do i=1,nclev call cpseti('PAI',i) call cpsetr('CLV',clev(i)) if (clev(i).ge.0.) then call cpseti('CLD',65535) else call cpsetc('CLD -- contour line dash pattern', + '$$''''$$''''') endif if (icolor.gt.0) call cpseti('CLC',0) enddo call cpsetc('ILT',' ') call cprect(f,nx,nx,ny,rwrk,lrwrk,iwrk,liwrk) if (icolor.le.0) then call cpcldr(f,rwrk,iwrk) call cplbdr(f,rwrk,iwrk) else call arinam(iama,liama) call cpclam(f,rwrk,iwrk,iama) call arscam(iama,xcra,ycra,lcra,iarea,igrp,10,colram) call cplbam(f,rwrk,iwrk,iama) call cplbdr(f,rwrk,iwrk) call cpcldm(f,rwrk,iwrk,iama,drawcl) endif call cpgetr('CIU',ciu) ! contour interval used if (pltmin.ge.pltmax) then call cpgetr('ZMN',zmin) ! field min call cpgetr('ZMX',zmax) ! field max else zmin = pltmin zmax = pltmax endif ! ! Add axes and labels: ! if (iht.le.0) then call labrect(xx,nx,yy,ny,'LATITUDE (DEG)','ZP (LN(P0/P))', + pltchsize) else call labrect(xx,nx,yy,ny,'LATITUDE (DEG)','HEIGHT (KM)', + pltchsize) endif ! ! Add extra right-hand axes in height and/or pmb: ! iyaxright = iyax-iyax/10*10 if (iyaxright.gt.0) + call yaxzpht(zp,fht,zp,nx,nzp,dzp,yy,ny,p0,iht,iyaxright, + vp,spval) end subroutine pltstrm end module plt