!
      module chapman_module
      use params_module,only: nlonp4,nlat,dz
      implicit none
      real,dimension(nlonp4,nlat) :: 
     |  chi,         ! solar zenith angle
     |  sin_chi,     ! sin(chi)
     |  cos_chi,     ! cos(chi)
     |  rt_sinchi,   ! sqrt(sin(chi))
     |  slt          ! local time
!
      contains
!-----------------------------------------------------------------------
      subroutine chapman(z,tn,ch4,h2,barm,vch4,vh2,vn2,scch4,
     |  sch2,scn2,lev0,lev1,lon0,lon1,lat)
      use params_module,only: dz
      use init_module,only: secs,sin_sundec,cos_sundec,idn,uthr
      use cons_module,only: pi,dphi,dlamda,re,rmass_ch4,rmass_h2,
     | rmass_n2
      use hist_module,only: modeltime ! for diag
!
! -----------------------------------------------------------------
! Titan TGCM modifications. Latest version: 9/21/04 S. W. BOUGHER
! Calculate line integrals for ch4,h2,n2:
!
! Titan TGCM modifications. Latest version: 4/26/05 S. W. BOUGHER
! Utilize Titan diurnal cycle (UT secs) for SLT calculation
! -----------------------------------------------------------------
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
!
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: 
     |  z ,   ! geopotential height (cm)
     |  tn,   ! neutral temperature (deg K)
     |  ch4,  ! methane (mmr)
     |  h2,   ! molecular hydrogen (mmr)
     |  barm  ! mean molecular weight
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) :: 
     |  vch4,  ! ch4 vertical integration
     |  vh2,   ! h2 vertical integration
     |  vn2,   ! n2 vertical integration
     |  scch4, ! ch4 slant column integration
     |  sch2,  ! h2 slant column integration
     |  scn2   ! n2 slant column integration
!
! Local:
      real :: rlat,coslat,sinlat,factor,rtpi,fmin,fmax
      integer :: k,i
      real :: rp(lev0:lev1,lon0:lon1), ! z+re
     |        ti(lev0:lev1,lon0:lon1)  ! tn at interfaces
      real :: 
     |        n2(lev0:lev1,lon0-2:lon1+2)  ! n2
      integer :: nlevs
      real,dimension(lev0:lev1,lon0:lon1) :: ! for diagnostic plotting
     |  chi_plt, slt_plt, idn_plt
!
! Set up sun related geometry according to current model date and time.
! Also set day/night index idn (init_module).
!
      nlevs = lev1-lev0+1 ! for addfsech
      rtpi = sqrt(pi)
      rlat = -.5*pi+(float(lat-1)+.5)*dphi
      coslat=cos(rlat)
      sinlat=sin(rlat)
      idn(lon0:lon1) = 1
!     Titan differences (09/17/04) ----------------------------------
!     -- slt and chi formulations same as Earth
!     -- slt calculated using 57600 Earth secs for 1-uthr for Titan 
!     -- chi limit is unique for Titan extended atmosphere  (off)
      do i=lon0,lon1
!       slt(i,lat)=amod(secs/3600.+(float(i-3)*dlamda+pi)*12./pi,24.)
        slt(i,lat)=amod(secs/57600.+(float(i-3)*dlamda+pi)*12./pi,24.)
        chi(i,lat)=acos(sin_sundec*sinlat+cos_sundec*coslat*
     |    cos(pi*(slt(i,lat)-12.)/12.))
        sin_chi(i,lat) = sin(chi(i,lat))
        cos_chi(i,lat) = cos(chi(i,lat))
        rt_sinchi(i,lat) = sqrt(sin_chi(i,lat))
!     ---- SZA>130 Extended Atmosphere of Titan
!       if (chi(i,lat) > 2.268919) idn(i) = 0
!     ---- SZA>105
!       if (chi(i,lat) > 1.8326)   idn(i) = 0
!     ---------------------------------------------------------------

        chi_plt(:,i) = chi(i,lat)
        slt_plt(:,i) = slt(i,lat)
        idn_plt(:,i) = idn(i)
      enddo                

!     write(6,"('chapman: lat=',i3,' modeltime=',4i4,
!    |  ' uthr=',f5.2,/,'slt=',(6e12.4))") lat,modeltime,uthr
!     write(6,"('chi=',(12f5.2))") chi(:,lat)
!     write(6,"('slt=',(12f5.2))") slt(:,lat)
!     write(6,"('idn=',(12i2))") idn(:)

!
!     call addfsech('CHI',' ',' ',chi_plt,lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('SLT',' ',' ',slt_plt,lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('IDN',' ',' ',idn_plt,lon0,lon1,nlevs,nlevs,lat)

!
! rp = z+re
      do k=lev0,lev1
        rp(k,:) = z(k,lon0:lon1)+re
      enddo
!     call addfsech('RP',' ',' ',rp,lon0,lon1,nlevs,nlevs,lat)
!
! ti is tn at interfaces:
      ti(1,:) = tn(lev1,lon0:lon1) ! tn bottom boundary is stored in top slot
      do k=lev0+1,lev1-1
        ti(k,:) = .5*(tn(k-1,lon0:lon1)+tn(k,lon0:lon1))
      enddo
      ti(lev1,:) = tn(lev1-1,lon0:lon1) ! nlevp1 <- nlev
!     call addfsech('TNI',' ',' ',ti,lon0,lon1,nlevs,nlevs,lat)
!
! Make n2 mmr:
      do k=lev0,lev1
        n2(k,lon0:lon1) = 1.-ch4(k,lon0:lon1)-h2(k,lon0:lon1)
      enddo
!     call addfsech('CHAP_CH4',' ',' ',ch4,lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('CHAP_H2',' ',' ',h2,lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('CHAP_N2',' ',' ',n2,lon0,lon1,nlevs,nlevs,lat)
!
! Calculate line integrals (vo[x] and sco[x] are output):
!
      call line_integ(ch4,rmass_ch4,vch4,scch4,ti,barm,
     |  cos_chi(lon0:lon1,lat),sin_chi(lon0:lon1,lat),
     |  rt_sinchi(lon0:lon1,lat),rp,idn(lon0:lon1),
     |  lon0,lon1,lev0,lev1,lat,1)
      call line_integ(h2,rmass_h2,vh2,sch2,ti,barm,
     |  cos_chi(lon0:lon1,lat),sin_chi(lon0:lon1,lat),
     |  rt_sinchi(lon0:lon1,lat),rp,idn(lon0:lon1),
     |  lon0,lon1,lev0,lev1,lat,2)
      call line_integ(n2,rmass_n2,vn2,scn2,ti,barm,
     |  cos_chi(lon0:lon1,lat),sin_chi(lon0:lon1,lat),
     |  rt_sinchi(lon0:lon1,lat),rp,idn(lon0:lon1),
     |  lon0,lon1,lev0,lev1,lat,3)
!
! Save column number densities:
!     call addfsech('VCH4',' ',' ',vch4(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('VH2',' ',' ',vh2(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('VN2',' ',' ',vn2(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!
! Save slant line integrals:
!     call addfsech('SCCH4',' ',' ',scch4(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('SCH2',' ',' ',sch2(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('SCN2',' ',' ',scn2(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!
      end subroutine chapman
!-----------------------------------------------------------------------
      subroutine line_integ(f,fmass,v,s,ti,barm,cos_chi,sin_chi,
     |  rt_sinchi,rp,idn,lon0,lon1,lev0,lev1,lat,if)
      use params_module,only: dz
      use cons_module,only: expz,grav,p0,avo,gask,pi,re
      implicit none
!
! Args:
      integer,intent(in) :: lon0,lon1,lev0,lev1,lat,if
      real,intent(in) :: fmass
      real,dimension(lev0:lev1,lon0:lon1),intent(in) :: 
     |  ti,        ! tn at interfaces
     |  rp         ! z+re         (S9)
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  f,         ! input species density (mmr)
     |  barm       ! mean molecular mass
      real,dimension(lon0:lon1),intent(in) ::
     |  cos_chi,   ! cosine solar zenith angle
     |  sin_chi,   ! sine solar zenith angle
     |  rt_sinchi  ! sqrt(sin_chi)
      integer,intent(in) :: 
     |  idn(lon0:lon1) ! day-night index
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  v,         ! output column density
     |  s          ! output slant-line integral
!
! Local:
      real,parameter :: big=1.e80
      integer :: k,i,nlevs
      real :: factor,rtpi,exparg
      real :: 
     |  rtrp(lev0:lev1,lon0:lon1),      ! sqrt(rp/2hp) (S6)
     |  yp  (lev0:lev1,lon0:lon1),      ! yp and ip    (S5)
     |  r2ig(lev0:lev1,lon0:lon1)       ! 2.*ig        (S4)
      real,external :: expo
!
      nlevs = lev1-lev0+1 ! for addfsech
!
! Top:
      factor = avo*p0*expz(lev1-1)*exp(-.5*dz)/(fmass**2*grav)
      do i=lon0,lon1
        v(:,i) = f(:,i)
      enddo
      where(v < 0.) v = 0.
      do i=lon0,lon1
        v(lev1,i) = factor*.5*(v(lev1-1,i)+v(lev1,i))*
     |                    barm(lev1,i)
      enddo
!
! Integrate from top down:
      do i=lon0,lon1
        do k=lev1-1,lev0,-1
          factor = avo*p0*expz(k)/(fmass*grav)*dz
          v(k,i) = v(k+1,i)+factor*v(k,i)
        enddo
      enddo
!
! Set up for slant-line integral:
      factor = fmass*grav/(2.*gask)
      rtpi = sqrt(pi)
      do i=lon0,lon1
        do k=lev0,lev1
          rtrp(k,i) = sqrt(rp(k,i)*factor/ti(k,i)) ! S6
          yp(k,i) = rtrp(k,i)*abs(cos_chi(i))      ! S5
          if (yp(k,i) >= 8.) then
            yp(k,i) = v(k,i)*rtpi*rtrp(k,i)*
     |        (0.56498823/(0.06651874+yp(k,i)))
          else
            yp(k,i) = v(k,i)*rtpi*rtrp(k,i)*
     |        ((1.0606963+0.5564383*yp(k,i))/((yp(k,i)+1.7245609)*
     |          yp(k,i)+1.0619896))
          endif
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
      factor = grav*fmass/gask
      do i=lon0,lon1
        do k=lev0,lev1
!
! Avoid exceeding max arg to exp on non-unicos platform:
          exparg = rp(k,i)*(1.-sin_chi(i))*factor/ti(k,i)

          if (idn(i)==1.and.exparg < 650.) then ! daytime
            r2ig(k,i) = 2.*v(k,i)*exp(exparg)*
     |        rtpi*rt_sinchi(i)*rtrp(k,i)
          else
            r2ig(k,i) = big
          endif

!         r2ig(k,i) = 2.*v(k,i)*expo(exparg,0)*
!    |      rtpi*rt_sinchi(i)*rtrp(k,i)
        enddo
      enddo
!
! Slant line integrals (0 if obscured by earth):
! Slant line integrals (0 if obscured by titan):
      do i=lon0,lon1
        do k=lev0,lev1
          if (cos_chi(i) >= 0.) then
            s(k,i) = big
          else
            s(k,i) = r2ig(k,i)-yp(k,i)
          endif
          if (rp(k,i)*sin_chi(i)-re < 0.) s(k,i) = big
          if (cos_chi(i) >= 0.) s(k,i) = yp(k,i)
        enddo
      enddo
      end subroutine line_integ
      end module chapman_module
