c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/frstz2.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine frstz2(z,thenot,fsx,f1,f2,ncp1)
c     ******************************************************************
c     subroutine frstz2
c
c     purpose
c     this routine is used to calculate the zero order source
c     functions for a parallel plane atmos. model.
c
c     method
c     the zero order source functions are computed by calculating
c     the intensity of the incident solar radiation at each source
c     function level using attenuation factors calculated for a
c     parallel plane atmosphere.
c
c     calling sequence
c        call frstz2(z,thnot,fsx,f1,f2,ncp1)
c
c        variable    type  i/o   description
c        --------    ----  ---   -----------
c
c        z(202)      r*8     o    zero order reduced source functions
c         thnot       r*8i     current solar zenith angle
c        fsx         r*8     o    direct flux reaching ground
c        f1          r*8     o    constant
c        f2          r*8     o    computational constant
c        ncp1        i*4     i    # levels in model atmos.
c
c     analysis and programming
c     k. f. klenk , p. m. smith sasc aug 77
c
c      modifications (date name purpose)
c      last modified by zia ahmad
c      purpose: to include the effect of molecular anisotropy
c
c     last modified 03/14/95...dave flittner
c     purpose: set pressure scale height used in gravity correction
c     to rayleigh scattering od.  Create new variable pscaleforg and
c     pass in common block consts.
c***********************************************************************
c
c
c
      implicit none
      include 'parameter.inc'
c     real *8 h(487),ps(487),xs(487),dxs(487),zs(487),cofx(4,487)
      integer ncp1
      real*8 z(202),thenot,fsx,f1,f2
      include 'consts.cmn'
      include 'thkns.cmn'
      include 'energy.cmn'
      include 'out.cmn'
c     new statement
      include 'depolt.cmn'
c     end of new statement
c local
      real*8 amuo,amuosq
      integer j
c
      amuo = dcos(thenot*cnvrt)
      if (thenot .eq. 90.0) amuo = 0.
      amuosq = amuo**2
c     new statements
      if(ipol.eq.0)then
          f1 = 0.1875d0*(1.d0 + amuosq)
          f2 = 0.1875d0*sq2*(1.d0 - amuosq)
      else
          f1=0.25d0*q1*(1.0d0+amuosq+2.0d0*q)
          f2=0.25d0*q1*delp*(1.0d0-amuosq)
      endif
c*****parallel plane atmosphere
         do 200 j=1,ncp1
         z(j)=dexp(-tt(j)/amuo)
200      continue
c
         fsx=amuo*z(ncp1)
      fs=fsx
c
c compute the energy input
c
      engyin=amuo
      return
      end