c
c------------------------------------------------------------------
c Begin file /home/sting/foster/libplt/contour.f
c------------------------------------------------------------------
c
      subroutine contour(zz,idimx,nx,ny,cint,cmin,cmax)
      dimension zz(idimx,ny)
      include 'cpwrk.h'
      data icls/13/
      external drawcl
c
c If cint=0, then let conpack pick contour levels, otherwise
c   force cint as contour level:
c
      if (cint.le.0.) then
        call cpseti('CLS -- contour level selector',icls)
        call cpsetr('CIS -- contour interval specifier',cint)
        call cpsetr('CMN -- contour minimum',cmin)
        call cpsetr('CMX -- contour maximum',cmax)
        call cpsetc('ILT - info label text',
     +    'MIN $ZMN$, MAX $ZMX$, INTERVAL $CIU$')
c
c If cint > 0, then hardwire min, max, and interval:
c   (cint > 0 from getinp, cmin < cmax, icls not used)
c
      else
        call cpseti('CLS -- contour level selector',1)
        call cpsetr('CIS -- contour interval specifier',cint)
        call cpsetr('CMN -- contour minimum',cmin)
        call cpsetr('CMX -- contour maximum',cmax)
c       write(6,"('contour hardwiring cmin=',e12.4,' cmax=',e12.4,
c    +    ' cint=',e12.4)") cmin,cmax,cint
      endif
c
c Initialize conpack and pick contour levels:
c (ncl = number of contour levels)
c
      call cprect(zz,idimx,nx,ny,rwrk,lrwrk,iwrk,liwrk)
      call cppkcl(zz,rwrk,iwrk)
      call cpgeti('NCL -- number of contour levels',ncl)
c
c Loop through contour levels:
c (set dashed pattern for levels < 0)
c
      do i=1,ncl
        call cpseti('PAI -- parameter array index',i)
        call cpgetr('CLV -- contour level',clev)
        if (clev.lt.0.)
     +    call cpsetc('CLD -- contour line dash pattern',
     +        '$$''''$$''''')
        call cpseti('CLU -- contour level use flag',3)
      enddo
      call arinam(iama,liama)
      call cplbam(zz,rwrk,iwrk,iama)
      call cpcldm(zz,rwrk,iwrk,iama,drawcl)
      call cplbdr(zz,rwrk,iwrk)
c
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine drawcl(xcs,ycs,ncs,iai,iag,nai)
      dimension xcs(*),ycs(*),iai(*),iag(*)
      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

c
c------------------------------------------------------------------
c Begin file /home/sting/foster/libplt/mkproj.f
c------------------------------------------------------------------
c
      subroutine mkproj(proj,viewport,cenlat,cenlon,rot,perimlat,eradii)
c
c Make an ezmap projection (type of projection in proj)
c
c On input:
c   proj = 'CE'   -> make cylindrical equidistant projection
c   proj = 'ST'   -> make stereographic polar projection
c   proj = 'SV'   -> make sat view projection
c   proj = 'MO'   -> make mollweide projection
c   viewport(4) = fractional coords left,right,bottom,top of
c                 desired location of projection
c   perimlat = latitude of perimeter for stereographic proj
c   utc = current ut (used to determine rotation for stereographic)
c
      character*(*) proj
      dimension viewport(4)
      data ncalls/0/
      save ncalls
c
      ncalls = ncalls+1
      if (ncalls.eq.1) then
        call mapstc('OU','CO')
        call mapsti('DO',1)
        call mapsti('G1',0)
        call mapsti('G2',2)
        call mapsti('VS',0)
      endif
c
c Cylindrical equidistant:
c
      if (proj.eq.'CE'.or.proj.eq.'MO') then
        call mapsti('EL',0)
        if (proj.eq.'MO') call mapsti('EL',1)
        call mappos(viewport(1),viewport(2),viewport(3),viewport(4))
        call mapset('MA',0.,0.,0.,0.)
        call mapint
        call maproj(proj,cenlat,cenlon,rot)
        call mapint
c
c Stereographic (polar):
c
      elseif (proj.eq.'ST') then
        call mapsti('EL',1)
        call mappos(viewport(1),viewport(2),viewport(3),viewport(4))
        call maproj(proj,cenlat,cenlon,rot)
        ang = 90.-abs(perimlat)
        call mapset('AN',ang,ang,ang,ang)
        call mapint
c
c Satellite view projection:
c
      elseif (proj.eq.'SV') then
        call mapsti('EL',1)
        call mapstc('OU','CO')
        call mapsti('DO',1)
        call mappos(viewport(1),viewport(2),viewport(3),viewport(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',0.1) ! distance (deg) between pts of grid lines
        call mapint
        call maproj(proj,cenlat,cenlon,rot)
        call mapint
c
c Bad proj string:
c
      else
        write(6,"('>>> mkproj: unsupported projection: ',a)") proj
      endif
      return
      end