!
      module set_vert
!
! Set up various vertical scaling for any plots requiring vertical
! coords on y-axis. Note these routines must be in a module, or
! have an explicit interface block in the calling routine, presumably
! because unallocated pointers are passed. 
!
      contains
!-------------------------------------------------------------------
      subroutine setvert(zpspec,gcmlev,nlev,dlev,nzprange,k0zp,k1zp,
     +  zprange,htspec,nhtscale,htscale,spv,iprint)
!
! Set up vertical scale(s) in zp and/or ht for plotting.
! (this may be called from any plot setup that requires an axis in 
!  vertical coordinates)
!
! On input:
!   zpspec(2) specifies bottom and top of requested zp range (if not spv)
!   htspec(3) specifies bottom, top and delta km of requested height 
!     scale (if not spv)
!   gcmlev(nlev) is array of valid zp
!   dlev is delta zp of gcmlev.
!   zprange and htscale are pointers to be allocated
!   spv = special value
!   iprint = print flag
!
! On output:
!   zprange(nzprange) is array of requested zp's (nzprange=0 if
!     zp was not requested)
!   if nzprange > 0, k0zp,k1zp are indices in gcmlev of zprange(1) and
!     zprange(nzprange)
!   htscale(nhtscale) is array of requested ht's (nhtscale=0 if
!     ht was not requested)
!
! Args:
      integer,intent(in) :: nlev,iprint
      real,intent(in) :: gcmlev(nlev),dlev,htspec(3),spv
      real,intent(inout) :: zpspec(2)
      integer,intent(out) :: k0zp,k1zp,nhtscale,nzprange
      real,pointer :: zprange(:),htscale(:)
!
! Locals:
      integer :: ier,k
!
! Externals:
      integer,external :: ixfind
!
! Return pointer zprange (zprange(nzprange)) according to zpspec:
! (zpspec = bottom, top (within gcmlev))
!
      nzprange = 1
      if (zpspec(1)/=spv.and.zpspec(2)/=spv) then
        if (zpspec(1) < gcmlev(1)) then
          if (iprint > 0) 
     +      write(6,"('Adjusted bottom of zp range from ',f8.3,
     +        ' to bottom of model at zp = ',f8.3)") zpspec(1),gcmlev(1)
          zpspec(1) = gcmlev(1)
        endif
        if (zpspec(2) > gcmlev(nlev)) then
          if (iprint > 0) 
     +      write(6,"('Adjusted top of zp range from ',f8.3,
     +        ' to top of model at zp = ',f8.3)")
     +      zpspec(2),gcmlev(nlev)
          zpspec(2) = gcmlev(nlev)
        endif
        if (zpspec(1) >= zpspec(2)) then
          write(6,"('>>> setvert: bad zprange=',2f9.3)") zpspec
          nzprange = 0
        endif
        if (nzprange > 0) then
          k0zp = ixfind(gcmlev,nlev,zpspec(1),dlev)
          k1zp = ixfind(gcmlev,nlev,zpspec(2),dlev)
          nzprange = k1zp-k0zp+1
          if (associated(zprange)) deallocate(zprange)
          allocate(zprange(nzprange),stat=ier)
          if (ier /= 0) 
     +      call allocerr(ier,"setvert allocating zprange")
          do k=k0zp,k1zp
            zprange(k-k0zp+1) = gcmlev(k)
          enddo
        endif
      else
        nzprange = 0
      endif
!
! Return pointer phtscale (htscale(nhtscale)) according to htspec:
! (htspec(3) = bottom, top, delta (kmx))
!
      nhtscale = 0
      if (htspec(1)/=spv.and.htspec(2)/=spv.and.
     +    htspec(3)/=spv) then
        if (htspec(2) <= htspec(1)) then
          if (iprint > 0) 
     +      write(6,"('>>> setvert: bad htspec (bottom must be ',
     +        '< top) htspec=',3f10.2)") htspec
          nhtscale = 0
        elseif (htspec(3) > htspec(2)-htspec(1)) then
          if (iprint > 0)
     +      write(6,"('>>> setvert: bad delta km for htscale = ',
     +        f10.2,' (must be < top-bot)')") htspec(3)
          nhtscale = 0
        else
          nhtscale = ifix((htspec(2)-htspec(1)) / htspec(3) + 1.0000001)
          if (associated(htscale)) deallocate(htscale)
          allocate(htscale(nhtscale),stat=ier)
          if (ier /= 0) 
     +      call allocerr(ier,"setvert allocating htscale")
          do k=1,nhtscale
            htscale(k) = htspec(1) + (k-1)*htspec(3)
          enddo
        endif
      endif
      end subroutine setvert
!-------------------------------------------------------------------
      subroutine sethscale(f,fhscale,nfhscale,htscale,nhtscale,
     +  hscale,nhscale)
      use fields,only: field
      implicit none
!
! Calculate a height scale to be used for a field which is calculated
!   only at a fixed height scale, i.e., oh-v and oh-b fields.
! On input:
!   f is a defined field type which was calculated only at its fixed
!     height range (f%nlev, f%dlev).
!   fhscale and hscale are pointers (output)
!   nhscale and nhscale are integers (output)
!   htscale(nhtscale) is the user requested height scale (input)
! On output:
!   fhscale(nfhscale) is the fixed height scale at which f was
!     calculated (set using f%lev, f%nlev,f%dlev)
!   hscale(nhscale) is the height scale to use when plotting the
!     field. This will use the bottom and top of the user's
!     requested range, but at the resolution of the field (overriding
!     the users requested resolution). Hscale may be all or partially
!     outside fhscale, so when data to be plotted is defined, you must
!     check for heights in hscale which are outside fhscale, and set 
!     data to spval accordingly.
!
! Args:
      type(field),intent(in) :: f
      real,pointer :: fhscale(:),hscale(:)
      integer,intent(out) :: nfhscale,nhscale
      integer,intent(in) :: nhtscale
      real,intent(in) :: htscale(nhtscale)
!
! Externals:
      real,external :: rnd
!
! Locals:
      integer :: ier,k
      real :: hbot,htop
!
! fhscale(nfhscale) is height scale allowable for current field:
!
      if (associated(fhscale)) deallocate(fhscale)
      nfhscale = f%nlev
      allocate(fhscale(nfhscale),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,"sethscale allocating fhscale") 
      do k=1,nfhscale
        fhscale(k) = f%lev(1)+float(k-1)*f%dlev
      enddo
!
! hscale(nhscale) is height scale to be used for current frame (assuming
!   field is height-only) using bottom and top of user requested htscale,
!   and vertical resolution of the field.
!
      if (associated(hscale)) deallocate(hscale)
      hbot = rnd(htscale(1),f%dlev/2.)
      htop = rnd(htscale(nhtscale),f%dlev/2.)
      nhscale = ifix((htop-hbot)/f%dlev+1.0000001)
      allocate(hscale(nhscale),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,"sethscale allocating hscale") 
      hscale(1) = hbot
      do k=2,nhscale
        hscale(k) = hscale(k-1)+f%dlev
      enddo
      return
      end subroutine sethscale
      end module set_vert