c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/getgasprf.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine getgasprf(igas)
c***********************************************************************
c
cccc
c     subroutine getgasprf
c
c     by: dave flittner
c
c     purpose-
c
c        Load other absorbing gas profiles.
c        Then calc the cumlative from top of atmo, pressure=1.0e-05 atm
c        to each altitude.  Then convert to pressure scale with 1 atm at 
c        the surface.
c
c        1. read the input profile, height vs concentration. Units are
c           height in km and concentration in #/cubic cm.
c        2. spline concentration to 0.25 km levels from surface (z=0 to z=81 km)
c        3. convert the altitude scale of column to pressure scale with the
c           2 arrays hhold and pshold.
c        4. spline column(pcolumn) at the 101 pressure levels
c
c     input/output variables -
c
c        variables      type   i/o  description
c        ---------      ----   ---  -----------
c     formal parameters
c        igas           i*4     i   index of the additional gases
c     part of common 'cgas2'
c        fgasprf(mgas) char*8   i   file name of the input profile
c        alphagas(mgas) r*8     i   gaseous absorption coefficient
c                                   Units=[cm*cm]
c        totgas(101,mgas) r*8   o   Total integrated gas igas from
c                                   top of atmosphere down to level
c
c last modified: March 13, 1998...def
c purpose: Change message about total column amount.
c
c last mod. Oct. 16, 1998...def
c purpose: use the new spline and splint that have mode switch
c
c last mod. Apr. 13,1999...def
c purpose: fix extrapolation of the column.
c
c     last mod. May 3, 1999...def
c     purpose: add option to compute addition absorbing gas profile using the
c     internal, implicit neutral density profile and a constant mixing ratio.
c     Currently it is set to do only o2 or o4 absorption and is controled by
c     the lo2abs and lo4abs logicals read in the ENV file.
c
c     last mod. May 6, 1999...def
c     purpose: set yp1 and ypn to use natrual spline.
c
c     last mod. Nov. 2, 1999...def
c     purpose: use .not. instead of .eqv. in logical testing
c
cccc
c***********************************************************************
c
      implicit none
      include 'parameter.inc'
      include 'consts.cmn'
      include 'cgas2.cmn'
      include 'switches.cmn'
      include 'atm_101.cmn'
c passed
      integer igas
c
c local
c
      integer n,mn,np,mp,i,j,iunt,mode,itype
      parameter (mn=82,mp=81*4+1)
      real*8 z(mn),con(mn),y2(mp),column(mp),zcolumn(mp),pcolumn(mp)
      real*8 yp1,ypn,tmp1,tmp2
      logical prtf,internalprf
      character*14 fname
c
c don't print files
c
      prtf=.FALSE.
      call get_lun(iunt)
c
      if(lo2abs.and.lo4abs.and.(igas.ge.ngas-2))then
        internalprf=.TRUE.
        if(igas.eq.ngas-2)then
          itype=1
        else if(igas.eq.ngas-1)then
          itype=2
        else
          itype=0
        endif
      else if(lo2abs.and.(.not.lo4abs).and.(igas.eq.ngas-1))then
        itype=1
        internalprf=.TRUE.
      else if(lo4abs.and.(.not.lo2abs).and.(igas.eq.ngas-1))then
        itype=2
        internalprf=.TRUE.
      else
        internalprf=.FALSE.
      endif
      if(internalprf)then
c
c compute oxygen number density profile using the pshold and hhold arrays
c
        call get_neutral_prf(hhold,pshold,mn,con,xo2,itype)
c let zcolumn = hhold
        np=mn
        do i=1,np
           zcolumn(i)=hhold(np-i+1)
           z(i)=hhold(i)
           pcolumn(i)=log(pshold(np-i+1))
        enddo
      else
c
c open gas profile and read
c
        open(iunt,file=fgasprf(igas),status='old')
        n=0
        do while (n.lt.mn)
           read(iunt,*,end=20)tmp1,tmp2
           if(tmp2.lt.0.0d0)then
             write(6,*)'getgasprf. negative concentraion read!'
             write(6,*)'getgasprf. for gas ',igas+1,' at alt ',tmp1
             write(6,*)'getgasprf. will stop now.'
             stop
           endif
           n=n+1
           z(n)=tmp1
           con(n)=tmp2
        enddo
20      close(iunt)
        if(n.eq.0)then
          write(6,*)'getgasprf. error reading profile for gas',igas+1
          write(6,*)'getgasprf. File name =',fgasprf(igas)
          write(6,*)'getgasprf. Will stop now'
          stop
        endif
c
        np=n
c let zcolumn = z
        do i=1,np
           zcolumn(i)=z(np-i+1)
        enddo
      endif
      if(internalprf.and.(itype.gt.0))then
        do i=1,np
           column(i)=con(np-i+1)
        enddo
      else
c compute the total column down to zcolumn(1) assuming an exponential dist
c convert from #/cubic cm to #/cm/cm/km
        j=np-1
        i=1
        if(con(j).ne.con(j+1))then
          column(i)=(con(j)-con(j+1))*
     &1.0d05*(z(j+1)-z(j))/log(con(j)/con(j+1))
        else
          column(i)=con(j+1)*1.0d05
        endif
c now integrate to the rest of the levels using a spline
        do i=2,np
           j=np-i+1
           if(con(j).ne.con(j+1))then
             column(i)=column(i-1)+(con(j)-con(j+1))*
     &1.0d05*(z(j+1)-z(j))/log(con(j)/con(j+1))
           else
             column(i)=column(i-1)+con(j)*1.0d05*(z(j+1)-z(j))
           endif
        enddo
c
c now convert the altitude scale of column to pressure scale with the
c 2 arrays hhold and pshold.
c now spline pshold(hhold)
        mode=2
        yp1=1.0e30
        ypn=1.0e30
        call spline(hhold,pshold,mn,yp1,ypn,y2,mode)
        do j=1,np
           if(zcolumn(j).gt.hhold(mn))then
             mode=6
             call spline(hhold,pshold,mn,yp1,ypn,y2,mode)
             call splint(hhold,pshold,y2,mn,zcolumn(j),pcolumn(j),
     &mode)
             mode=2
             call spline(hhold,pshold,mn,yp1,ypn,y2,mode)
           else 
             call splint(hhold,pshold,y2,mn,zcolumn(j),pcolumn(j),mode)
           endif
           pcolumn(j)=log(pcolumn(j))
        enddo
      endif
      if(prtf)then
        if(igas+1.lt.10)then
          write(fname(1:11),'(a6,i1,a4)')'column',igas+1,'.out'
          open(iunt,file=fname(1:11),status='unknown')
        else
          write(fname(1:12),'(a6,i2,a4)')'column',igas+1,'.out'
          open(iunt,file=fname(1:12),status='unknown')
        endif
        do j=1,np
           write(iunt,'(1p3e12.4)')zcolumn(j),
     &exp(pcolumn(j)),column(j)
        enddo
        close(iunt)
        if(igas+1.lt.10)then
          write(fname,'(a8,i1,a4)')'othergas',igas+1,'.out'
          open(iunt,file=fname(1:13),status='unknown')
        else
          write(fname,'(a8,i2,a4)')'othergas',igas+1,'.out'
          open(iunt,file=fname(1:14),status='unknown')
        endif
      endif
c
c now spline column(pcolumn) at the 101 pressure levels
c
      mode=2
      call spline(pcolumn,column,np,yp1,ypn,y2,mode)
      do i = 1, 101
         if(p_101lg(i).lt.pcolumn(1))then
           if((column(2).gt.column(1)))then
c use exponetial to extrapolate
             totgas(i,igas)=log(column(1))+
     &((exp(p_101lg(i))-exp(pcolumn(1)))/
     &(exp(pcolumn(2))-exp(pcolumn(1)))*
     &(log(column(2))-log(column(1))))
           totgas(i,igas)=exp(totgas(i,igas))
          else
           totgas(i,igas)=column(1)
          endif
         else
           call splint(pcolumn,column,y2,np,p_101lg(i),
     &totgas(i,igas),mode)
         endif
         if(prtf)write(iunt,*)exp(p_101lg(i)),totgas(i,igas)
      enddo
      if(prtf)close(iunt)
c
c subtract off the amount above pmin, so that totgas(1)=0.0
c
      do i=101,1,-1
         totgas(i,igas)=totgas(i,igas)-totgas(1,igas)
      enddo
c     write(6,1000)igas+1,totgas(101,igas),totgas(101,igas)/Nlosch
1000  format(' getgas. Total column amount for gas ',i2,'=',1pe10.3,
     &'[#/cm^2] or',1pe10.3,'[atm-cm]')
1010  format(' getgas. Total column amount for gas ',i2,'=',1pe10.3,
     &'[#/cm^2] or',1pe10.3,'[atm-cm]',/,
     &       '         Absorption vertical optical depth=',1pe10.3)
      return
      end
c***********************************************************************