c
      subroutine yaxzpht(hts1d,hts2d,gcmzp,id1,kmx,dzp,yy,ny,p0,
     +  iht,iyax)
c
c Add extra non-linear y-axis on right in height or pressure:
c
      real yy(ny),hts1d(kmx),hts2d(id1,kmx),gcmzp(kmx)	! input
      real yaxht(kmx),pmbht(ny),yaxmb(kmx)		! local
c
c For pressure (mb) right hand y-axis: note strings pfcode (with
c function codes) must match values of presslab (see yaxright.f)
c
      parameter(npresslab=14,nhtlab=6)
      dimension presslab(npresslab),htlab(nhtlab)
      character*8 pfcode(npresslab)
      data presslab
     +  / 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/
      data pfcode  
     +  /'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'/
c
      call pcsetc('FC','$')
c
c If zp is on left y-axis (iht <= 0), then add non-linear ht y-axis
c
      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     write(6,"('yaxzpht: p0=',e12.4,' kmx=',i2,' nzp=',i2,
c    +  ' yaxmb(1-nzp)=',/(6(1pe12.4)))") p0,kmx,nzp,(yaxmb(k),k=1,nzp)
c
c Offsets to new right y-axis from existing (unlabeled) right y-axis:
c
      chsize = 0.	! default will be .014
      offtoax = 0.	! default will be .015
      offtolab = 0.	! default will be offtoax+.09
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
        if (iyax.gt.1) then
          chsize = .012
          offtolab = offtoax+.085
        endif
        call yaxright(yy,yaxht,ny,htlab,nhtlab,'(I3)',1,pfcode,
     +    nminor,'AVE HEIGHT (KM)',offtoax,offtolab,chsize)
c
c Optionally add pressure(mb) axis in addition to height:
c
        if (iyax.gt.1) then
          offtolab = 0.
          offtoax = .12
          nminor = 0
          call yaxright(yy,yaxmb,ny,presslab,npresslab,' ',0,pfcode,
     +      nminor,'PRESSURE (MB)',offtoax,offtolab,chsize)
        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)
      endif
      return
      end