subroutine firstz(h,ps,z,thenot,fracin,lmax,ncp1)
c
c***********************************************************************
c
cccc
c     subroutine firstz
c
c     purpose-
c        correct the optical depths calculated for the layers of a plane
c        parallel atmosphere for the sphericity of the actual atmosphere
c
c        method-
c      1.using the optical depths of the layers of the standard atmos.
c        calculate the slant optical path of the solar rays thru them
c        using the chapman function.
c      2.then using spline interpolation calculate the slant optical
c        paths in the model atmosphere.
c
c        calling sequence-
c        call firstz(h,ps,z,thenot,fracin,fs,f1,f2,lmax,ncp1)
c
c      variables     type   i/o description
c      ---------     ----   --- -----------
c      h(487)        r*8    i   height of layers of standard atmosphere
c                               from earth center (units of r)
c      ps(487)        r*8   i   rayleigh optical thickness of
c                               each layer of std. atmosphere
c      z(202)         r8    o   total optical path of solar ray to
c                               reach each layer of the model atmos.
c      thenot         r*8       solar zenith angle
c      fracin         r*8   i   layer*r (radius of earth in units of
c                               layer thickness)
c      fs             r*8   o   cos(thenot)*z(202)
c      f1             r*8   o   3/16(cos(thenot)**2)
c      f2             r*8   o   3*sqrt(2)/16(sin(thenot)**2)
c      lmax           i*4   i   # layers of standard atmos.
c      ncp1           i*4   i   # layers in model atmos.
c
c     external references
c        omerf
c        splset
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.
cccc
c**********************************************************************
c
c     ******************************************************************
c
c        using the slant path optical thicknesses and chapman function
c        this subroutine calculates dave's z matrices for primary
c        scattering.
c
c     ******************************************************************
      implicit integer*4(i-n),real*8 (a-h,o-z)
      real *8 h(487),ps(487),zs(487),cofx(4,487)
      real*8 z(202)
c
      include "consts.inc"
      include "out.inc"
      include "thkns.inc"
      include "chpmn.inc"
c
c*****calculate trig functions to be used in do loops
c
      amuo = dcos(thenot*cnvrt)
      if (thenot .eq. 90.0) amuo = 0.
      amuosq = amuo**2
      sn = dsqrt(1.-amuosq)
      f1 = 0.1875*(1. + amuosq)
      f2 = 0.1875*sq2*(1. - amuosq)
c
c*****calculate slant opt. path of solar rays to reach
c*****layer(i) of standard atmosphere
c
      do 200 i = 1, lmax
      zs(i) = 0.
      raycon = sn*h(i)
      if (i .eq. 1) go to 190
      duma= dsqrt((h(1) + raycon)*(h(1) - raycon))
c
c*****add slant optical paths thru each layer(j)
c*****lying between layer(1) and layer(i)
c
      do 180 j = 2, i
      dumb= dsqrt((h(j) + raycon)*(h(j) - raycon))
      path = fracin*(duma - dumb)
      zs(i) = zs(i) + path*dxs(j)
  180 duma = dumb
c
c*****find the slant opt. path to reach layer 1 from the top of
c*****the atmosphere- using chapman function approximation
c
  190 amu = raycon/h(1)
      amu = dsqrt(1. - amu**2)
      chpp = omerf(sqchp*amu, amu, chpn)
      chxx = omerf(sqchx*amu, amu, chxn)
c
c*****add the previous results to obtain total slant path zs to reach
c*****each layer i.
c
      zs(i) = zs(i) + chpp*ps(1) + chxx*(dxs(1) - ps(1))
  200 continue
c
c*****fit a spline between xs=log(vert. path) and zs=log(slant path)
c*****through the standard atmosphere
c
      do 215 i = 1, lmax
  215 zs(i) = dlog(zs(i))
      call splset(xs,zs,cofx,lmax)
c
c*****for each layer of model atmosphere find two layers in the
c*****standard atmosphere straddling it and find total slant optical
c*****path (z) to reach each layer of model atmosphere-
c*****   using spline interpolation
c
      z(1) = 1.0
      j = 2
      lmaxm1 = lmax - 1
      do 230 i = 1, lmaxm1
      if (ttl(j) .le. xs(i) .or. ttl(j) .gt. xs(i+1)) go to 230
  220 dum1 = xs(i+1) - ttl(j)
      dum2 = ttl(j) - xs(i)
      dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)*
     1dum2**2 + cofx(4,i))
      dum3 = dexp(dum3)
      z(j) = dexp(-dum3)
      j = j + 1
      if (j .gt. ncp1) go to 240
      if (ttl(j) .gt. xs(i) .and. ttl(j) .le. xs(i+1)) go to 220
  230 continue
c
c*****find z for the bottom layer by extrapolation if it was not
c*****obtained during the above interpolation
c
  240 if (ttl(ncp1) .le. xs(lmax)) go to 245
      dum1 = xs(lmax) - ttl(ncp1)
      dum2 = ttl(ncp1) - xs(lmaxm1)
      i = lmaxm1
      dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)*
     1dum2**2 + cofx(4,i))
      dum3 = dexp(dum3)
      z(ncp1) = dexp(-dum3)
  245 continue
      fs = amuo*z(ncp1)
c     write (6,'(i4,e12.4)') (i,z(i),i=1,ncp1)
      return
      end