c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/cmp_gc.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine cmp_gc()
c ----------------------------------------------------------------
      implicit none
c
c compute the gravity correction factor for the molecular optical depth
c
      include 'parameter.inc'
      include 'consts.cmn'
      include 'atm_101.cmn'
      include 'switches.cmn'
c local
      integer i
      real*8 Hp,lnp0
      real*8 introvrzprsqr
      external introvrzprsqr
c
c gc_type= 0 - no gravity correction
c          1 - explicitly old one (1-pscaleforg*p_101lg)**2
c          2 - new routine for entire atmo
c          3 - mimic old with new routine (lnp0=0,Hp=6.95,gc_factor(z=0)=1.0)
c          4 - partially mimic old with new routine (lnp0=0,Hp=6.95),
c              but don't impose gc_factor(z=0)=1.0.
c
c
      if(gc_type.eq.0)then
        do i=1,101
           gc_factor(i)=1.0d0
        enddo
      else if(gc_type.eq.1)then
        do i=1,101
           gc_factor(i)=(1.0d0-pscaleforg*p_101lg(i))**2
        enddo
      else if((gc_type.ge.2).and.(gc_type.le.4))then
        i=1
        Hp=(alt_101(1)-alt_101(2))/(p_101lg(2)-p_101lg(1))
        lnp0=p_101lg(1)+(alt_101(1)/Hp)
        if((gc_type.eq.3).or.(gc_type.eq.4))then
          lnp0=0.0
          Hp=6.95d0
        endif
        gc_factor(1)=p_101(1)*introvrzprsqr(lnp0,p_101lg(1),r,Hp)
        if(jprint(6).ne.0)write(6,1000)'cmp_gc. gc_factor ',
     &alt_101(i),p_101(i),gc_factor(i)/p_101(i)
1000    format(a20,1P3e14.6)
        do i=2,101
           Hp=(alt_101(i-1)-alt_101(i))/(p_101lg(i)-p_101lg(i-1))
           lnp0=p_101lg(i)+(alt_101(i)/Hp)
           if((gc_type.eq.3).or.(gc_type.eq.4))then
              lnp0=0.0
              Hp=6.95d0
           endif
           gc_factor(i)=gc_factor(i-1)+
     &(p_101(i)*introvrzprsqr(lnp0,p_101lg(i),r,hp)-
     &p_101(i-1)*introvrzprsqr(lnp0,p_101lg(i-1),r,hp))
           if(jprint(6).ne.0)write(6,1000)'cmp_gc. gc_factor ',
     &alt_101(i),p_101(i),gc_factor(i)/p_101(i)
        enddo
        do i=1,101
           gc_factor(i)=gc_factor(i)/p_101(i)
        enddo
        if(gc_type.eq.3)then
c normalize such that at the surface the correction factor is 1.
          do i=1,101
             gc_factor(i)=gc_factor(i)/gc_factor(101)
          enddo
        endif
      endif
      return
      end