c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/firstzmuop.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine firstzmuop(h,ps,z,thenot,fracin,lmax,ncp1,theta,phi)
c
c***********************************************************************
c
cccc
c     subroutine firstzmuop
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        to the line of sight specified by theta and phi.
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 firstzmuop(h,ps,z,thenot,fracin,fs,lmax,ncp1,theta,phi)
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      lmax           i*4   i   # layers of standard atmos.
c      ncp1           i*4   i   # layers in model atmos.
c      theta          r*8   i   view angle at the surface
c      phi            r*8   i   azimuthal angle at the surface
c
c     external references
c        omerf
c        splset
c
c    last modified 03/10/95...dave flittner
c    purpose: change the solar zenith angle along the line of sight
c    using the input solar zenith angle as the solar zenith angle
c    at the surface.  The direction of the line of sight will at the
c    surface and be taken from emu(1) (the first scan angle), and the 
c    azimuth will be taken also at the surface from the first azimuth angle.
c
c    last modified 06/27/95...dave flittner
c    purpose: modify for solar zenith angle greater than 90.0.  The 
c    logical lsubhoriz is set to true if solar zenith angle is greater
c    than 90.  Now will have points along the zenith that are in the
c    shadow (zs=0).  Because of this the log of zs can not be taken.
c    Must spline to zs.
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 none
      include 'parameter.inc'
c passed
      integer lmax,ncp1
      real *8 h(487),ps(487),z(202),theta,phi,thenot,fracin
c commons
      include 'chpmn.cmn'
      include 'consts.cmn'
      include 'out.cmn'
      include 'thkns.cmn'
c functions
      real*8 omerf
      external omerf
c local
      real*8 zs(487),cofx(4,487),gamma,raycon
      real*8 cosgamma,singamma,cosphi,costhenotp,sinthenotp
      real*8 amuo,amuosq,sn,sintheta,xshadow,dumb,duma,path,dum1
      real*8 dum2,dum3,amu,chpp,chxx
      integer ishadow,i,ij,j,lmaxm1
      logical lsubhoriz
c
      ishadow=0
c
c*****calculate trig functions to be used in do loops
c
      amuo = dcos(thenot*cnvrt)
      if (thenot .eq. 90.0) amuo = 0.0d0
      amuosq = amuo**2
      sn = dsqrt(1.0d0-amuosq)
c
c*****calculate slant opt. path of solar rays to reach
c*****layer(i) of standard atmosphere
c
      cosphi=dcos(phi*cnvrt)
      sintheta=dsin(theta*cnvrt)
      do 200 i = 1, lmax
c calc. the angle, gamma, between the local zenith and the ref. zenith
         gamma=theta-dasin(1.0d0/h(i)*sintheta)/cnvrt
         cosgamma=dcos(gamma*cnvrt)
         singamma=dsin(gamma*cnvrt)
         costhenotp=amuo*cosgamma-sn*singamma*cosphi
         sinthenotp=dsqrt(1.0d0-costhenotp*costhenotp) 
         zs(i) = 0.0d0
         raycon = sinthenotp*h(i)
c to keep the solar zenith angle the same along the line of sight,
c set raycon equal to sn*h(i).
c i.e. raycon = sn*h(i)
         if(costhenotp.le.0.0d0)then
           lsubhoriz=.TRUE.
           if(raycon.le.1.0d0)then
             do ij=i,lmax
                zs(ij)=0.0d0
             enddo
             xshadow=xs(i)
             ishadow=i
             goto 210
           endif
         else
           lsubhoriz=.FALSE.
         endif
         if (i .gt. 1)then
           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
           if(lsubhoriz)then
             j=i+1
             do while(h(j).gt.raycon)
                dumb= dsqrt((h(j) + raycon)*(h(j) - raycon))
c factor of 2 because there is a contribution from this layer at 2
c points along the line of the direct solar beam.
                path = 2.0d0*fracin*(duma - dumb)
                zs(i) = zs(i) + path*dxs(j)
                duma=dumb
                j=j+1
             enddo
c now the tangent layer
             path = 2.0d0*fracin*(duma)
             zs(i) = zs(i) + path*dxs(j)
           endif
         endif
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.0d0 - 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
  210 if(ishadow.gt.0)then
        do i = 1, ishadow-1
           zs(i) = dlog(zs(i))
        enddo
c extrapolate so can use the spline function since if set shadow points
c to zero, the spline will oscillate.  Will use the xshadow value
c to determine if the grid point is actually in the shade.
        do i=ishadow,lmax
           zs(i)=zs(ishadow-1)+(xs(i)-xs(ishadow-1))*
     &(zs(ishadow-1)-zs(ishadow-2))/(xs(ishadow-1)-xs(ishadow-2))
        enddo
      else
        do 215 i = 1, lmax
  215   zs(i) = dlog(zs(i))
      endif
      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.0d0
      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))
      if(ishadow.gt.0)then
        if(i+1.ge.ishadow)then
          do ij=j,ncp1
             z(j)=0.0d0
          enddo
          goto 245
        endif
      endif
      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))
      if(ishadow.gt.0)then
        if(lmax.ge.ishadow)then
          z(ncp1)=0.0d0
        else
          dum3 = dexp(dum3)
          z(ncp1) = dexp(-dum3)
        endif
      else
        dum3 = dexp(dum3)
        z(ncp1) = dexp(-dum3)
      endif
  245 continue
      fs = amuo*z(ncp1)
c	write(6,'(a16,1E14.8,2f6.1,2E16.9)')'firstzmuop. z',
c     &z(ncp1),theta,phi,dexp(-dexp(zs(i))),dexp(-dexp(zs(i+1)))
c     write (6,'(i4,e12.4)') (i,z(i),i=1,ncp1)
      return
      end