!-----------------------------------------------------------------------
      subroutine hnmf2(ht,ne,hmf2,nmf2,fof2,lev0,lev1,lon0,lon1,lat)
!
! Find the NmF2, foF2, and hmF2 based on the electron profile.
! btf: This is taken from an idl procedure provided by Liying.
! fof2 added April, 2012.
!
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,intent(in),dimension(lev0:lev1,lon0:lon1) :: ht,ne
      real,intent(out),dimension(lon0:lon1) :: hmf2,nmf2,fof2
!
! Local:
      integer :: k,kx,i
      real :: h(3),p(3),h12,h22,h32,deltx,atx,ax,btx,bx,ctx,cx
!
      do i=lon0,lon1
        kx = 0
        kloop: do k = lev1-1,2,-1
          if (ne(k,i) >= ne(k-1,i).and.ne(k,i) >= ne(k+1,i)) then
            kx = k
            exit kloop
          endif
        enddo kloop
        if (kx==0) then
          write(6,"('>>> hnmf2: lat=',i4,' could not find kx -- ',
     |      'setting hmf2=nmf2=fof2=0')") lat
          hmf2(:) = 0. ; nmf2(:) = 0. ; fof2(:) = 0.
          return
        endif
        h = (/ht(kx-1,i),ht(kx,i),ht(kx+1,i)/)
        p = (/ne(kx-1,i),ne(kx,i),ne(kx+1,i)/)
        h12 = h(1)*h(1)
        h22 = h(2)*h(2)
        h32 = h(3)*h(3)
        deltx=h12*h(2)+h22*h(3)+h32*h(1)-h32*h(2)-h12*h(3)-h22*h(1)
        atx=p(1)*h(2)+p(2)*h(3)+p(3)*h(1)-h(2)*p(3)-h(3)*p(1)-h(1)*p(2)
!
! Periodic points are not set in ht and ne, so they are zero
! at i=1,2. Setting hmf2and nmf2 zero here is ok, since they
! are saved on histories at i=3,nlon. Without this check, code 
! will stop here w/ divide by zero deltx or zero atx when debug 
! flags are set.
!
        if (deltx == 0..or.atx == 0.) then
          hmf2(:) = 0. ; nmf2(:) = 0. ; fof2(:) = 0.
          cycle
        endif
        ax=atx/deltx
        btx=h12*p(2)+h22*p(3)+h32*p(1)-h32*p(2)-h12*p(3)-h22*p(1)
        bx=btx/deltx
        ctx=h12*h(2)*p(3)+h22*h(3)*p(1)+h32*h(1)*p(2)-h32*h(2)*p(1)-
     |    h12*h(3)*p(2)-h22*h(1)*p(3)
        cx=ctx/deltx
        hmf2(i)=-(bx/(2.*ax))
        nmf2(i)=-((bx*bx-4.*ax*cx)/(4.*ax))
        fof2(i)=sqrt(nmf2(i)*1.e6/1.24e10)
      enddo ! i=lon0,lon1
!     write(6,"('hnmf2: lat=',i4,' hmf2 min,max=',2e12.4,
!    |  ' nmf2 min,max=',2e12.4,' fof2 min,max=',2e12.4)") 
!    |  lat,minval(hmf2),maxval(hmf2),minval(nmf2),maxval(nmf2),
!    |  minval(fof2),maxval(fof2)
      end subroutine hnmf2
!-----------------------------------------------------------------------
!
! This is Liying's original IDL procedure hnmf2.pro. This was converted
! to fortran in subroutine hnmf2 above.
!
! PRO hnmf2, X, Y, XOUX, YOUX
! ;
! ; Find the NmF2 and hmF2 based on the electron profile
! ;
! ; input
! ;     X: height
! ;     Y: Ne profile
! ;
! ; output
! ;     XOUX: hmF2
! ;     YOUX: NmF2
! ;
! nz=n_elements(x)
! for k=nz-2,1,-1 do begin
!   if ((y[k] ge y[k-1]) and (y[k] ge y[k+1])) then break
! endfor
! h=[x[k-1],x[k],x[k+1]]
! p=[y[k-1],y[k],y[k+1]]
! 
! h12=h(0)*h(0)
! h22=h(1)*h(1)
! h32=h(2)*h(2)
! DELTX=h12*h(1)+h22*h(2)+h32*h(0)-h32*h(1)-h12*h(2)-h22*h(0)
! ATX=p(0)*h(1)+p(1)*h(2)+p(2)*h(0)-h(1)*p(2)-h(2)*p(0)-h(0)*p(1)
! AX=ATX/DELTX
! BTX=h12*p(1)+h22*p(2)+h32*p(0)-h32*p(1)-h12*p(2)-h22*p(0)
! BX=BTX/DELTX
! CTX=h12*h(1)*p(2)+h22*h(2)*p(0)+h32*h(0)*p(1)-h32*h(1)*p(0)-h12*h(2)*p(1)-h22*h(0)*p(2)
! CX=CTX/DELTX
! XOUX=-(BX/(2.*AX))
! YOUX=-((BX*BX-4.*AX*CX)/(4.*AX))
!
! RETURN
! END
!-----------------------------------------------------------------------

