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
      implicit integer*4(i-n),real*8 (a-h,o-z)
c     real *8 h(487),ps(487),xs(487),dxs(487),zs(487),cofx(4,487)
      real *8 z(202)
c
      include "consts.inc"
      include "thkns.inc"
      include "chpmn.inc"
      include "out.inc"
c
c     new statement
      include "depolt.inc"
c     end of new statement
c
      amuo = dcos(thenot*cnvrt)
      if (thenot .eq. 90.0) amuo = 0.
      amuosq = amuo**2
      sn = dsqrt(1.-amuosq)
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
      return
      end