c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/get_neutral_prf.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
        subroutine get_neutral_prf(hhold,pshold,nhold,
     &n_neutral,x_neutral,itype)
c
c
c compute the neutral density profile given the contant mixingratio
c
c input:
c hhold - altitude profile
c phhold - pressure profile
c nhold - # of points in the profile
c x_neutral - mixing ratio
c itype - type of profile output using the mixing ratio x_neutral
c         0 - concentration
c         1 - column amount of concentration to the pressure level
c         2 - column amount of squared concentration to the pressure level
c output;
c n_neutral - number density of the profile
c
c p = g * m * int(n(z)dz)
c p= [atm] * 101325 [N/m^2]/[atm] 
c g= 9.81 [m/s/s]
c m= 28.964e-03/6.022e23 [kg/particle]
c int(n(z)dz) * 1.0e09 = [particle/cm/cm] if n(z) = [particle/cm^3] &
c                                            z in km.
c
        implicit none
        include 'consts.cmn'
        include 'switches.cmn'
c passed
        integer nhold,itype
        real*8 hhold(nhold),pshold(nhold),n_neutral(nhold),x_neutral
c local
        integer i
        real*8 dp,Hp_top,gcorrect,factor,n_bottom,n_top
c
c constant for convertion of pressure [atm] to number density [1/cc]
c        factor=101325/(9.81*28.964e-03/6.022e23*1.0e09)
        factor=p0/(g0*mair/Na*1.0d09)
        gcorrect=1.0d0
c for the top layer assume that the pressure scale height and density 
c are the same
        if(gc_type.gt.0)gcorrect=(1.0d0+hhold(nhold)/r)**2
        Hp_top=(hhold(nhold)-hhold(nhold-1))/
     &log(pshold(nhold-1)/pshold(nhold))
        n_bottom=pshold(nhold)*factor/Hp_top*gcorrect*x_neutral
        if(itype.eq.0)then
          n_neutral(nhold)=n_bottom
        else if(itype.eq.1)then
          n_neutral(nhold)=pshold(nhold)*factor*gcorrect*
     &x_neutral*1.0d05
        else if(itype.eq.2)then
          n_neutral(nhold)=Hp_top/2.0d0*n_bottom*n_bottom*1.0d05
        endif 
        i=nhold
c       write(6,1000)i,hhold(i),pshold(i),n_bottom,n_neutral(i)
        n_top=n_bottom
        do i=nhold-1,1,-1
           Hp_top=(hhold(i+1)-hhold(i))/log(pshold(i)/pshold(i+1))
           if(gc_type.gt.0)
     &gcorrect=(1.0d0+0.5*(hhold(i)+hhold(i+1))/r)**2
           dp=(pshold(i)-pshold(i+1))
           n_bottom=n_top+dp*factor/Hp_top*gcorrect*x_neutral
           if(itype.eq.0)then
             n_neutral(i)=n_bottom
           else if(itype.eq.1)then
             n_neutral(i)=n_neutral(i+1)+dp*factor*gcorrect*
     &x_neutral*1.0d05
           else if(itype.eq.2)then
             n_neutral(i)=n_neutral(i+1)+
     &Hp_top/2.0d0*(n_bottom*n_bottom-n_top*n_top)*1.0d05
           endif 
c          write(6,1000)i,hhold(i),pshold(i),n_bottom,n_neutral(i)
           n_top=n_bottom
        enddo
c
1000    format(i4,f6.2,1P3E12.4)
        end