!
      subroutine getoh(f,nf,flat,fo2,fo1,fn2,fh,fo3,fho2,foh,
     +  h,imx,kmx,nlat,lat,ixlat,ixlon)
!
! Define OH vibrational levels and/or OH emission bands at
!   current latitude lat:
! If lat<0, then get oh fields (vertical columns) only at given 
!   lat,lon indices ixlat,ixlon
! On input, flat data is defined at latitude index abs(lat)
!
      use fields,only: field,nohalt,oh_alt,set_ohalt,nohvlev,nohband
      use hist,only: history
      use proc,only: spval
      use input,only: ibohv_watts
      implicit none
!
! Args:
      integer,intent(in) :: nf,imx,kmx,nlat,lat,ixlat,ixlon
      type(history),intent(in) :: h
      real,intent(in) :: flat(imx,kmx,h%nflds)
      real,intent(in) :: fo2(imx,kmx),fo1(imx,kmx),fn2(imx,kmx),
     +  fh(imx,kmx),fo3(imx,kmx),fho2(imx,kmx),foh(imx,kmx)
      type(field),intent(inout) :: f(nf)
!
! Locals:
      integer :: i,ix,ixt,ixz,ier,len,lo,hi,ilat
      integer,parameter :: ideltav=2
      real :: rho(kmx)
      real :: tn_ohv(nohalt),   xo2_ohv(nohalt),  xo_ohv(nohalt), 
     +        xn2_ohv(nohalt),  xh_ohv(nohalt),   xo3_ohv(nohalt), 
     +        xho2_ohv(nohalt), xoh_ohv(nohalt),  rho_ohv(nohalt)
      real :: fohv(nohvlev,nohalt)	! output of ohrad (cm3)
      real :: fohb(60,nohalt)		! output of ohrad (kilo Rayleighs)
      integer :: ixfohv(nf), ifohv(nf)
      integer :: ixfohb(nf), ifohb(nf)
!
! Externals:
      integer,external :: ixfindc,ixohband
!
! ilat is index to f()%data(:,ilat,:,:)
! If lat < 0, then we want field only at latitude index ixlat,
!   so skip out now if lat < 0 and not at desired latitude this call:
!
      ilat = abs(lat)
      if (lat < 0) then	! caller wants data only at ixlat,ixlon
        if (ilat /= ixlat) return
      endif
!
! OH band emissions are dependent on t, o2, o, n2, h, o3, ho2, oh:
      ixz  = ixfindc(h%fnames,h%nflds,'Z       ')
      ixt  = ixfindc(h%fnames,h%nflds,'TN      ')
!
! Define indices to needed ohv and ohb fields, and corresponding 
!   indices to fohv,fohb ouputs from ohrad:
!
      ixfohv = 0
      ifohv = 0
      ixfohb = 0
      ifohb = 0
      fields_loop: do ix=1,nf
        if (.not.f(ix)%derived) cycle fields_loop
        len = len_trim(f(ix)%fname8)
        if (trim(f(ix)%type) == 'OH-VIB') then
          ixfohv(ix) = ix
          read(f(ix)%fname8(len:len),"(i1)") ifohv(ix)
        elseif (trim(f(ix)%type) == 'OH-BAND') then
          ixfohb(ix) = ix
          ifohb(ix) = ixohband(f(ix)%fname8,hi,lo)
          if (ibohv_watts > 0) f(ix)%units = "WATTS/CM3/SEC"
        endif
      enddo fields_loop
!
! Interpolate dependencies to oh height scale:
!     subroutine intloc(gcmin,gcmht,kmx,hts,nhts,loght,
!    +  gcmout,idim1,ndim1,idim2,ier,spval,iprnt)
!
      lonloop: do i=1,imx
        if (lat < 0) then	! only at ixlat,ixlon
           if (i /= ixlon) cycle lonloop
        endif
        call intloc(flat(i,:,ixt),flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,tn_ohv,1,1,nohalt,ier,spval,0)
        call intloc(fo2(i,:),flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,xo2_ohv,1,1,nohalt,ier,spval,0)
        call intloc(fo1(i,:),flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,xo_ohv,1,1,nohalt,ier,spval,0)
        call intloc(fn2(i,:),flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,xn2_ohv,1,1,nohalt,ier,spval,0)
        call intloc(fh(i,:),flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,xh_ohv,1,1,nohalt,ier,spval,0)
        call intloc(fo3(i,:),flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,xo3_ohv,1,1,nohalt,ier,spval,0)
        call intloc(fho2(i,:),flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,xho2_ohv,1,1,nohalt,ier,spval,0)
        call intloc(foh(i,:),flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,xoh_ohv,1,1,nohalt,ier,spval,0)
        rho(:) = fo2(i,:)+fo1(i,:)+fn2(i,:)
        call intloc(rho,flat(i,:,ixz),kmx,oh_alt,nohalt,
     +    0,rho_ohv,1,1,nohalt,ier,spval,0)
!
!      subroutine ohrad(kmbot,kmtop,mvl,nkm,TA,SO2,SO,SN2,SH,SO3,SHO2,
!    +   SOH0,SM,ideltav,dohv,bohv,iwatts,spval,iprint)
!
        call ohrad(ifix(oh_alt(1)),ifix(oh_alt(nohalt)),nohvlev,nohalt,
     +    tn_ohv,xo2_ohv,xo_ohv,xn2_ohv,xh_ohv,xo3_ohv,xho2_ohv,xoh_ohv,
     +    rho_ohv,ideltav,fohv,fohb,ibohv_watts,spval,0)
        do ix=1,nf
          if (ixfohv(ix) > 0)
     +      f(ixfohv(ix))%data(i,ilat,:) = fohv(ifohv(ix),:)  
          if (ixfohb(ix) > 0) then
            f(ixfohb(ix))%data(i,ilat,:) = fohb(ifohb(ix),:)  
          endif
        enddo
      enddo lonloop	! i=1,imx
      return
      end