!
      subroutine getna(f,nf,flat,h,imx,kmx,lat)
!
! Calculate Na species requested in f(nf) for current latitude lat:
!
! The following 13 (mxnaf) fields are partitioned from NAT by napart.f:
!   (see local na_names(mxnaf))
!     'NaS     ','NaO     ','NaO3    ','NaO2    ','NaOH    ',
!     'NaCO3   ','NaHCO3  ','NaS+    ','NaN2+   ','NaCO2+  ',
!     'NaH2O+  ','NaO+    ','NaEMIS  '
! These have the following 13 dependencies:
!   TN, O2, H2O, CO2, H, H2, O1, O3, NE, O2+ [NO+], O+, N2 NAT 
!   (NAT is "total sodium" from the history)
!   (since NO+ is not available on histories, napart will use 
!    NO+ = (Ne - O+ - O2+))
!
      use fields,only: field,mxnaf
      use hist,only: history
      use proc,only: gcmlat,gcmlon,nlon,dlon,spval
      implicit none
!
! Args:
      integer,intent(in) :: nf,imx,kmx,lat
      type(field),intent(inout) :: f(nf)
      type(history),intent(in) :: h
      real,intent(in) :: flat(imx,kmx,h%nflds)
!
! Locals:
      integer :: i,ix,nna,ndep,ier,iday
      real :: slt,fn2(kmx)
      character(len=8) :: reqnames(mxnaf) ! names of requested sp
      character(len=8),allocatable :: 
     +  depnames(:)			! names of Na dependencies
      integer :: ixna(mxnaf)		! indices to requested na fields
      integer,allocatable :: ixdep(:)	! indices to dependences in flat
      real,allocatable :: fna(:,:)	! output from napart
      real :: fmin,fmax
!
! Externals:
      integer,external :: ixfindc
      real,external :: fslt
!
      reqnames = ' '	! names of requested Na fields
      nna = 0		! number of requested Na fields
      ndep = 0		! number of dependencies (same for all Na sp)
      do i=1,nf
        if (f(i)%fname8(1:2)=='Na') then
          nna = nna+1
          reqnames(nna) = f(i)%fname8
          ixna(nna) = i
          if (nna==1) then	! assume all Na sp have same dependencies
            if (associated(f(i)%fneed)) then
              ndep = size(f(i)%fneed)
              allocate(depnames(ndep),stat=ier)
              if (ier /= 0) call allocerr(ier,
     +          'getna allocating dependency names for Na')
              depnames = f(i)%fneed	! dependency names
              allocate(ixdep(ndep),stat=ier)
              if (ier /= 0) call allocerr(ier,
     +          'getna allocating dependency indices for Na')
            endif
          endif
        endif
      enddo
      if (nna == 0) then
        write(6,"('>>> getna: no Na species requested -- returning')")
        return
      endif
      if (ndep == 0) then
        write(6,"('>>> getna: no dependencies?? ndep=',i2)") ndep
        return
      endif
!
! Find indices to dependencies in flat: 
! (assume all Na sp have same dependencies). 
! These must be in this order for the call to napart (see fset_known):
!   TN, O2, H2O, CO2, H, H2, O1, O3, NE, O2+ O+, N2 NAT 
!
      do i=1,ndep
        ixdep(i) = ixfindc(h%fnames,h%nflds,depnames(i))
!
! OP is alias for O+ on netcdf histories:
        if (i==11.and.ixdep(i)==0) 
     |    ixdep(i) = ixfindc(h%fnames,h%nflds,'OP      ')
!
! O2P is alias for O2+ on netcdf histories:
        if (i==10.and.ixdep(i)==0) 
     |    ixdep(i) = ixfindc(h%fnames,h%nflds,'O2P     ')
        if (ixdep(i)==0.and.depnames(i)/='N2      ') then
          write(6,"('>>> WARNING getdep: cannot find Na dependency '
     +      ,a,' in history fields.')") depnames(i)
        endif
      enddo
      deallocate(depnames)
!
! Output of napart is fna(kmx,nna)
      allocate(fna(kmx,nna),stat=ier) 
      if (ier /= 0) call allocerr(ier,
     +  'getna allocating fna data for Na species.')
      iday = h%iyd-h%iyd/1000*1000
!
! Lon loop:
      do i=1,imx
        slt = fslt(slt,h%ut,gcmlon(i),1)
        fn2(:) = max(.00001,(1.-flat(i,:,ixdep(2))-flat(i,:,ixdep(7))))
!
!     subroutine napart(tn,xo2,xh2o,xco2,xh,xh2,xo1,xo3,xne,
!    +  xo2p,xop,xn2,xnat,ncol,iden,zsb,zst,npart,names,fout,
!    +  iday,slt,glat,glon,spval,idebug)
! (note flat densities are mmr, so pass napart iden=0)
!
        call napart(flat(i,:,ixdep(1)),		! tn
     +              flat(i,:,ixdep(2)),		! o2
     +              flat(i,:,ixdep(3)),		! ho2
     +              flat(i,:,ixdep(4)),		! co2
     +              flat(i,:,ixdep(5)),		! h
     +              flat(i,:,ixdep(6)),		! h2
     +              flat(i,:,ixdep(7)),		! o1
     +              flat(i,:,ixdep(8)),		! o3
     +              flat(i,:,ixdep(9)),		! ne
     +              flat(i,:,ixdep(10)),	! o2+
     +              flat(i,:,ixdep(11)),	! o+
     +              fn2,			! n2
     +              flat(i,:,ixdep(13)),	! nat
     +    kmx,0,h%zpb,h%zpt,nna,reqnames,fna,iday,slt,gcmlat(lat),
     +    gcmlon(i),spval,0)
!
! Transfer fields to f():
        do ix = 1,nna
          f(ixna(ix))%data(i,lat,:) = fna(:,ix)
        enddo
      enddo
      deallocate(ixdep)
      deallocate(fna)
      end subroutine getna