!
! Make various map projections, contour fields and plot vectors:
!
      subroutine pltce(map,vp,spv)
!
! Make cylindrical equidistant map projection w/ contoured field:
!
      use plt
      use mk_maps,only: mapproj
      use proc,only: dlat,gcmlat,nlat
      use input,only: multiplt,ipltrowcol
      implicit none
!
! Args:
      type(mapproj),intent(in) :: map
      real,intent(in) :: vp(4),spv
!
! Locals:
      integer,parameter :: nlonlab=13, nlatlab=7
      integer :: latlab(nlatlab)=(/-90,-60,-30,0,30,60,90/)
      integer :: lonlab(nlonlab),i,iclr,incx=3,incy=3
      real :: xmid,dum,chsize,yboffset,dlatlab
      character*8 vectype
      integer :: jlat0,jlat1,ny
!
! Externals:
      real,external :: convlon
      integer,external :: ixfind
!
! Set y-axis to plot according to map%pltlat0,1:
! pltlat0,1 are from namelist read fmap_latminmax(2)
! map%yy, map%ny, and map%data are currently dimensioned
!   and defined at the full gcmlat(nlat) range.
!
      jlat0 = ixfind(map%yy,map%ny,map%pltlat0,dlat)
      jlat1 = ixfind(map%yy,map%ny,map%pltlat1,dlat)
      if (jlat0<0.or.jlat0>map%ny.or.jlat1<0.or.jlat1>map%ny.or.
     |  jlat0>=jlat1) then
        write(6,"('>>> pltmaps: bad jlat0,1=',2i4,
     |    ' map%pltlat0,1=',2f9.3)") jlat0,jlat1,map%pltlat0,
     |    map%pltlat1
        write(6,"('>>> will resort to full latitude range..')")
        jlat0 = 1
        jlat1 = map%ny
      endif
      ny = jlat1-jlat0+1
      if (map%pltlat0 /= gcmlat(1).or.map%pltlat1 /= gcmlat(nlat)) then
        latlab(1) = map%pltlat0
        latlab(nlatlab) = map%pltlat1
        dlatlab = (map%pltlat1-map%pltlat0)/6.
        do i = 2,6
          latlab(i) = latlab(i-1)+dlatlab
        enddo
!       write(6,"('pltce: pltlat0,1=',2f8.2,' dlatlab=',f8.2,
!    |    ' gcmlat0,1=',2f8.2)") map%pltlat0,map%pltlat1,dlatlab,
!    |    gcmlat(jlat0),gcmlat(jlat1)
!       write(6,"('latlab=',7i4)") latlab
      endif

!     write(6,"('pltce: jlat0,1=',2i4,' ny=',i4)") 
!    |  jlat0,jlat1,ny
!
      xmid = 0.5*(vp(1)+vp(2))
      call cpseti('SET',0)
      call cpseti('MAP',1)
      call cpsetr('XC1',map%xx(1))
      call cpsetr('XCM',map%xx(map%nx))
      call cpsetr('YC1',map%yy(jlat0))
      call cpsetr('YCN',map%yy(jlat1))
!
! Make projection:
! subroutine mkproj(proj,vp,cenlat,cenlon,rot,perimlat,eradii,icont)
      call mkproj(map%proj,vp,0.,map%cenlon,0.,dum,dum,0,
     |  map%pltlat0,map%pltlat1)
!
! Contour field:
! If map%vectype=='     ' (blank),    -> contour only (no vectors)
! If map%vectype=='+UNVN' or '+UIVI', -> add vectors to contours
! If map%vectype=='UNVN ' or 'UIVI ', -> vectors only (no contour)
! (if vectors only, use vec sum (map%data) for zmin,max)
!
      if (len_trim(map%vectype)==0.or.map%vectype(1:1)=='+') then
        call contour(map%data(:,jlat0:jlat1),map%nx,map%nx,ny)
      else
        call fminmax(map%data(:,jlat0:jlat1),
     |    size(map%data(:,jlat0:jlat1)),zmin,zmax)
      endif
!
! Add vector arrows:
!     subroutine pltvec(u,v,id1,id2,x0,x1,y0,y1,incx,incy,
!    +  yboffset,iclr,iveclab,veclab,spv)
!
      if (len_trim(map%vectype) > 0) then
!       yboffset = .75	! fraction of viewport height
        yboffset = .80	! fraction of viewport height
        iclr = iwhite
        if (icolor>0.and.index(map%vectype,'+U')>0) iclr = iblack
        if (index(map%vectype,'UN')>0) vectype = 'UN+VN'
        if (index(map%vectype,'UI')>0) vectype = 'UI+VI'
        call pltvec(map%u,map%v,map%nx,ny,map%xx(1),map%xx(map%nx),
     +    map%yy(1),map%yy(ny),incx,incy,yboffset,iclr,labelv,
     +    vectype,spv,1)
      endif
!
! Add continental outlines:
      if (map%icontinents > 0) then
        iclr = iwhite
        if (icolor > 0.and.
     +    map%vectype(1:1)/='U'.and.map%vectype(1:1)/='V') 
     +    iclr = iblack
        call mpseti('C5 -- color index for continents',iclr)
        call maplot
      endif
!
! X,Y Axes:
!
      do i=1,nlonlab ! every 30 deg
        lonlab(i) = ifix(convlon(map%cenlon-180.+(i-1)*30.,180))
      enddo
      chsize = .017*(vp(2)-vp(1))
!     if (multiplt > 0 .and. ipltrowcol(2)==1) 
!    +  chsize = .030*(vp(4)-vp(3))
      call labaxes(nlonlab,lonlab,3,map%xlab,
     +             nlatlab,latlab,3,map%ylab,chsize)
      end subroutine pltce
!-----------------------------------------------------------------------
      subroutine pltst(map,vp,ut,glats,inc,nlat,spv)
!
! Make polar stereographic map projection w/ contoured field:
!
      use input,only: isattrack
      use plt
      use mk_maps,only: mapproj
      implicit none
!
! Args:
      type(mapproj),intent(in) :: map
      real,intent(in) :: vp(4),ut,spv
      integer,intent(in) :: nlat
      integer,intent(in) :: inc(nlat)
      real,intent(in) :: glats(nlat)
!
! Locals:
      real :: dlat,gperimlat,yc1,ycn,rot,plat,dum,xmid,fmin,fmax,
     +  yboffset
      integer :: ixperim,lat0,lat1,nlats,i,j,ier,iclr
      integer :: noinc(nlat)
      real,allocatable :: fpol(:,:),upol(:,:),vpol(:,:),vsum(:,:)
      character*8 vectype
!
! Externals:
      integer,external :: ixfind
!
! Prepare conpack:
      xmid = 0.5*(vp(1)+vp(2))
      call cpseti('SET',0)
      call cpseti('MAP',1)
      call cpsetr('XC1',map%xx(1))
      call cpsetr('XCM',map%xx(map%nx))
!
! Define lats, rot, etc for hemisphere and perimlat:
      dlat = glats(2)-glats(1)
      ixperim = ixfind(glats,nlat,map%perimlat,dlat)
      gperimlat = glats(ixperim)	! nearest grid perimlat
      if (gperimlat > 0. .and. gperimlat < map%perimlat) then
        gperimlat = gperimlat+dlat
        ixperim = ixperim+1
      endif
      if (gperimlat < 0. .and. gperimlat > map%perimlat) then
        gperimlat = gperimlat-dlat
        ixperim = ixperim-1
      endif
      if (gperimlat > 0.) then	! north polar
        yc1 = gperimlat
        ycn = glats(nlat)
        lat0 = ixperim
        lat1 = nlat
        rot = -ut*15.
        plat = 90.
      else			! south polar
        yc1 = glats(1)
        ycn = gperimlat
        lat0 = 1
        lat1 = ixperim
        rot = ut*15.-180.
        plat = -90.
      endif
      nlats = lat1-lat0+1
!
! Make stereographic projection:
      call cpsetr('YC1',yc1)
      call cpsetr('YCN',ycn)
      call mkproj('ST',vp,plat,0.,rot,map%perimlat,dum,0,
     |  map%pltlat0,map%pltlat1)
!
! Define plot array only at needed latitudes:
      allocate(fpol(map%nx,nlats),stat=ier)
      if (ier /= 0) call allocerr(ier,"pltst allocating fpol") 
      noinc = 0	! array
      call mkfpol(map%data,fpol,map%nx,map%ny,lat0,lat1,noinc,spv)
!
! Contour field:
! If map%vectype=='     ' (blank),    -> contour only (no vectors)
! If map%vectype=='+UNVN' or '+UIVI', -> add vectors to contours
! If map%vectype=='UNVN ' or 'UIVI ', -> vectors only (no contour)
!
      if (len_trim(map%vectype)==0.or.map%vectype(1:1)=='+')
     +  call contour(fpol,map%nx,map%nx,nlats)
!
! Add vector arrows:
!     subroutine pltvec(u,v,id1,id2,x0,x1,y0,y1,incx,incy,
!    +  yboffset,iclr,iveclab,veclab,spv)
!
      if (len_trim(map%vectype) > 0) then
        allocate(upol(map%nx,nlats),stat=ier) 
        if (ier /= 0) call allocerr(ier,"pltst allocating upol") 
        allocate(vpol(map%nx,nlats),stat=ier) 
        if (ier /= 0) call allocerr(ier,"pltst allocating vpol") 
        call mkfpol(map%u,upol,map%nx,map%ny,lat0,lat1,inc,spv)
        call mkfpol(map%v,vpol,map%nx,map%ny,lat0,lat1,inc,spv)
        yboffset = .20	! fraction of viewport height
        iclr = iwhite
        if (icolor>0.and.index(map%vectype,'+U')>0) iclr = iblack
        if (index(map%vectype,'UN')>0) vectype = 'UN+VN'
        if (index(map%vectype,'UI')>0) vectype = 'UI+VI'
        call pltvec(upol,vpol,map%nx,nlats,map%xx(1),map%xx(map%nx),
     +    map%yy(lat0),map%yy(lat1),1,1,yboffset,iclr,labelv,vectype,
     +    spv,1)
!
! If vectors only, use zmin,max of vector sum:
        if (trim(map%vectype)=='UNVN'.or.trim(map%vectype)=='UIVI')then
          allocate(vsum(map%nx,nlats),stat=ier)
          if (ier /= 0) call allocerr(ier,"pltst allocating vsum") 
          vsum = spv	! init array
          do j=1,nlats
            do i=1,map%nx
              if (upol(i,j)/=spv.and.vpol(i,j)/=spv) 
     +          vsum(i,j)=sqrt(upol(i,j)**2+vpol(i,j)**2)
            enddo
          enddo 
          call fminmax(vsum,size(vsum),zmin,zmax)
          deallocate(vsum)
        endif
        deallocate(upol) 
        deallocate(vpol) 
      endif
!
! Add continental outlines:
      if (map%icontinents > 0) then
        iclr = iwhite
        if (icolor > 0.and.
     +    map%vectype(1:1)/='U'.and.map%vectype(1:1)/='V') 
     +    iclr = iblack
        call mpseti('C5 -- color index for continents',iclr)
        call maplot
      endif
!
! Draw circle around perimlat and label w/ local time:
      call labpolar(vp(1),vp(2),vp(3),vp(4),.020,gperimlat)
!
! Release local allocatable:
      deallocate(fpol) 
      return
      contains
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Sub labpolar is internal to pltst:
      subroutine labpolar(pmapl,pmapr,pmapb,pmapt,charsz,perimlat)
      use mk_utlat,only: issltmap
      implicit none
!
! Put local time labels on a polar plot:
! (perimlat used only to determine hemisphere)
!
! Args:
      real,intent(in) :: pmapl,pmapr,pmapb,pmapt,charsz,perimlat
!
! Locals:
      character(len=2) :: sltlab
      integer :: i
      integer :: icent(12) = (/(-1,i=1,3),0,(1,i=1,5),0,(-1,i=1,2)/)
      integer :: mdlab=2
      real,parameter :: pi=3.14156
      real :: dgrid=30.,dtor=pi/180.
      real :: vl,vr,vb,vt,wl,wr,wb,wt,theta,cx,cy,r,x1,y1,x2,y2,
     +  szfs,tlenmjr,tlenmnr,xx,yy
      integer :: lty,mlab,minr,ng,ii
!
      call getset(vl,vr,vb,vt, wl,wr,wb,wt, lty)
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      theta = 0.
!
! (CX,CY) = center of circle
! R = radius of circle
      cx = pmapl + (pmapr-pmapl)*0.5
      cy = pmapb + (pmapt-pmapb)*0.5
      r = pmapr - cx
!
! Draw circle:
      do i=0,359
        x1 = cx + r*cos(float(i)*dtor)
        y1 = cy + r*sin(float(i)*dtor)
        x2 = cx + r*cos(float(i+1)*dtor)
        y2 = cy + r*sin(float(i+1)*dtor)
        call line(x1,y1,x2,y2)
      enddo
      if(issltmap > 0 .or. isattrack > 0) then
        call set(vl,vr,vb,vt, wl,wr,wb,wt, lty)
        return
      endif
      szfs = .015
      if (charsz.gt.0.) szfs = charsz
      szfs = szfs*(vr-vl)
      call plch(pmapl+0.5*(pmapr-pmapl),pmapb-.05*(vr-vl),
     +  'LOCAL TIME',szfs,0.,0.)
!
! Set up tics and numeric labels around the circle:
      mlab = 6
      if (perimlat.lt.0.) mlab = 18
      write(sltlab,"(i2)") mlab
      ng = 360./dgrid
      theta = 0.
      x1 = cx + r*cos(theta)
      y1 = cy + r*sin(theta)
      tlenmjr = .025*(pmapr-pmapl)
      tlenmnr = 0.5*tlenmjr
      minr = 1 
      do i=1,ng
        x2 = cx + (r-tlenmjr)*cos(theta)
        y2 = cy + (r-tlenmjr)*sin(theta)
        call line(x1,y1,x2,y2)
        xx = cx + (r+tlenmjr)*cos(theta)
        yy = cy + (r+tlenmjr)*sin(theta)
        if (i.eq.1) xx = xx-.018
        call plch(xx,yy,sltlab,szfs,0.,float(icent(i)))
        do ii=1,minr-1
          theta = theta + dgrid/float(minr) * dtor
          x1 = cx + r*cos(theta)
          y1 = cy + r*sin(theta)
          x2 = cx + (r-tlenmnr)*cos(theta)
          y2 = cy + (r-tlenmnr)*sin(theta)
          call line(x1,y1,x2,y2) 
        enddo
        theta = theta + dgrid/float(minr) * dtor
        if (perimlat.lt.0.) then
          if (mlab-mdlab.lt.0) mlab = mlab-mdlab+26
          if (mlab-mdlab.ge.0) mlab = mod(mlab-mdlab,24)
        endif
        if (perimlat.gt.0.) mlab = mod(mlab+mdlab,24)
        write(sltlab,"(i2)") mlab
        x1 = cx + r*cos(theta)
        y1 = cy + r*sin(theta)
      enddo
      call set(vl,vr,vb,vt, wl,wr,wb,wt, lty)
      return
      end subroutine labpolar
      end subroutine pltst
!-----------------------------------------------------------------------
      subroutine pltsv(map,vp,inc,nlat,spv)
!
! Make satellite view map projection w/ contoured field:
!
      use plt
      use mk_maps,only: mapproj
      implicit none
!
! Args:
      integer,intent(in) :: nlat,inc(nlat)
      type(mapproj),intent(in) :: map
      real,intent(in) :: vp(4),spv
!
! Locals:
      real :: xmid,dgrid=30.,dum,yboffset
      integer :: iclr,ier
      character*8 vectype
      real,allocatable :: upol(:,:),vpol(:,:)
!
! Prepare conpack:
      xmid = 0.5*(vp(1)+vp(2))
      call cpseti('SET',0)
      call cpseti('MAP',1)
      call cpsetr('XC1',map%xx(1))
      call cpsetr('XCM',map%xx(map%nx))
      call cpsetr('YC1',map%yy(1))
      call cpsetr('YCN',map%yy(map%ny))
!
! Set up map projection with center at cenlatlon or cenlatslt:
      call mapstr('GR',dgrid)
      call mkproj('SV',vp,map%cenlat,map%cenlon,0.,dum,map%eradii,0,
     |  map%pltlat0,map%pltlat1)
!
! Contour field:
! If map%vectype=='     ' (blank),    -> contour only (no vectors)
! If map%vectype=='+UNVN' or '+UIVI', -> add vectors to contours
! If map%vectype=='UNVN ' or 'UIVI ', -> vectors only (no contour)
!
      if (len_trim(map%vectype)==0.or.map%vectype(1:1)=='+') then
        call contour(map%data,map%nx,map%nx,map%ny)
      else
        call fminmax(map%data,size(map%data),zmin,zmax)
      endif
!
! Add vector arrows:
!     subroutine pltvec(u,v,id1,id2,x0,x1,y0,y1,incx,incy,
!    +  yboffset,iclr,iveclab,veclab,spv)
!
      if (len_trim(map%vectype) > 0) then
        allocate(upol(map%nx,map%ny),stat=ier) 
        if (ier /= 0) call allocerr(ier,"pltsv allocating upol") 
        allocate(vpol(map%nx,map%ny),stat=ier) 
        if (ier /= 0) call allocerr(ier,"pltsv allocating vpol") 
        call mkfpol(map%u,upol,map%nx,map%ny,1,map%ny,inc,spv)
        call mkfpol(map%v,vpol,map%nx,map%ny,1,map%ny,inc,spv)
        yboffset = .20	! fraction of viewport height
        iclr = iwhite
        if (icolor>0.and.index(map%vectype,'+U')>0) iclr = iblack
        if (index(map%vectype,'UN')>0) vectype = 'UN+VN'
        if (index(map%vectype,'UI')>0) vectype = 'UI+VI'
        call pltvec(upol,vpol,map%nx,map%ny,map%xx(1),map%xx(map%nx),
     +    map%yy(1),map%yy(map%ny),1,1,yboffset,iclr,labelv,vectype,
     +    spv,1)
        deallocate(upol) 
        deallocate(vpol) 
      endif
!
! Add continental outlines and grid:
      iclr = iwhite
      if (icolor > 0.and.
     +  map%vectype(1:1)/='U'.and.map%vectype(1:1)/='V') 
     +  iclr = iblack
      if (map%icontinents > 0) then
        call mpseti('C5 -- color index for continents',iclr)
        call maplot
      endif
      call mpseti('C2 -- color index for grid',iclr)
      call mpseti('C1 -- color index for perimeter',iclr)
      call mapgrd
      call maplbl
      return
      end subroutine pltsv
!-----------------------------------------------------------------------
      subroutine pltmo(map,vp,spv)
!
! Make mollweide map projection w/ contoured field:
!
      use plt
      use mk_maps,only: mapproj
      implicit none
!
! Args:
      type(mapproj),intent(in) :: map
      real,intent(in) :: vp(4),spv
!
! Locals:
      real :: xmid,dgrid=30.,dum,yboffset
      integer :: idashpat(16)=(/1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0/)
      integer :: iclr,incx=3,incy=3
      character*8 vectype
!
! Externals:
      integer,external :: ipat
!
! Prepare conpack:
      xmid = 0.5*(vp(1)+vp(2))
      call cpseti('SET',0)
      call cpseti('MAP',1)
      call cpsetr('XC1',map%xx(1))
      call cpsetr('XCM',map%xx(map%nx))
      call cpsetr('YC1',map%yy(1))
      call cpsetr('YCN',map%yy(map%ny))
!
! Set up map projection:
      call mkproj('MO',vp,map%cenlat,map%cenlon,0.,dum,dum,0,
     |  map%pltlat0,map%pltlat1)
!
! Contour field:
! If map%vectype=='     ' (blank),    -> contour only (no vectors)
! If map%vectype=='+UNVN' or '+UIVI', -> add vectors to contours
! If map%vectype=='UNVN ' or 'UIVI ', -> vectors only (no contour)
! (if vectors only, use vec sum (map%data) for zmin,max)
!
      if (len_trim(map%vectype)==0.or.map%vectype(1:1)=='+') then
        call contour(map%data,map%nx,map%nx,map%ny)
      else
        call fminmax(map%data,size(map%data),zmin,zmax)
      endif
!
! Add vector arrows:
!     subroutine pltvec(u,v,id1,id2,x0,x1,y0,y1,incx,incy,
!    +  yboffset,iclr,iveclab,veclab,spv)
!
      if (len_trim(map%vectype) > 0) then
        yboffset = .30	! fraction of viewport height
        iclr = iwhite
        if (icolor>0.and.index(map%vectype,'+U')>0) iclr = iblack
        if (index(map%vectype,'UN')>0) vectype = 'UN+VN'
        if (index(map%vectype,'UI')>0) vectype = 'UI+VI'
        call pltvec(map%u,map%v,map%nx,map%ny,map%xx(1),map%xx(map%nx),
     +    map%yy(1),map%yy(map%ny),incx,incy,yboffset,iclr,labelv,
     +    vectype,spv,1)
      endif
!
! Add continental outlines and grid:
      iclr = iwhite
      if (icolor > 0.and.
     +  map%vectype(1:1)/='U'.and.map%vectype(1:1)/='V') 
     +  iclr = iblack
      if (map%icontinents > 0) then
        call mpseti('C5 -- color index for continents',iclr)
        call maplot
      endif
      call mpseti('C1 -- color index for perimeter',iclr)
      call mpseti('C2 -- color index for grid',iclr)
      call mapstr('GR -- grid resolution',dgrid)
      call mapsti('DO -- dotted outlines',1)
      call mapsti('DA -- grid dashline',ipat(idashpat))
      call mapdrw
      return
      end subroutine pltmo
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine mkfpol(fglb,plt,nx,ny,lat0,lat1,inc,spv)
      implicit none
!        
! Transfer fglb(:,lat0-lat1) to plt(:,1-nlat), where nlat = lat1-lat0+1
! If inc(1) > 0, then inc is array of increments in longitude in which
!   to define plt, E.g., if inc(18)=8 then define plt at every 8th
!   longitude around the 18th latitude, values of plt at the other
!   longitudes get spv special value (and are therefore not plotted)
!   (inc is used to limit number of vectors plotted onto polar
!    stereographic projections)
!
! Args:
      integer,intent(in) :: nx,ny,lat0,lat1
      real,intent(in) :: fglb(nx,ny),spv
      real,intent(out) :: plt(nx,lat1-lat0+1)
      integer,intent(in) :: inc(ny)
!
! Locals:
      integer :: j,i
!
      if (inc(1) <= 0) then     ! do not use inc
        do j=lat0,lat1
          plt(:,j-lat0+1) = fglb(:,j)
        enddo
      else                      ! use inc
        do j=lat0,lat1
          do i=1,nx
            if (mod(i,inc(j)) == 0) then
              plt(i,j-lat0+1) = fglb(i,j)
            else
              plt(i,j-lat0+1) = spv
            endif
          enddo
        enddo
      endif
      return
      end subroutine mkfpol