!
      module chapman_module
      use params_module,only: nlonp4,nlat,dz
      use addfld_module,only: addfld
      implicit none
!
! Chi: Solar zenith angle (global). Used by chapman, heatnirco2, oplus, 
!   and settei. Calculated from local time (slt), and sun's declination
!   sundec, sin_sundec, cos_sundec (See init module -- these are updated 
!   daily by sub advance_day (advance.F)).
!
      real,dimension(nlonp4,nlat) ::
     |  chi          ! solar zenith angle
      integer,dimension(nlonp4,nlat) ::
     |  idn          ! day/night index
      real,parameter :: big=1.e80
!
      contains
!-----------------------------------------------------------------------
      subroutine chapman(z,tn,o2,o1,ox,n2,o3,co2,no,barm,vo2,vo1,vn2,
     |  vo3,vco2,vno,sco2,sco1,scn2,sco3,scco2,scno,lev0,lev1,lon0,lon1,
     |  lat)
!
! Calculate vertical and slant-line integrals.
! 8/4/04 btf: added NO.
!
      use params_module,only: dz
      use init_module,only: secs,sin_sundec,cos_sundec,uthr
      use cons_module,only: pi,dphi,dlamda,re,rmass_o2,rmass_o1,
     |  rmass_n2,rmass_o3,rmass_co2,rmass_no
      use hist_module,only: modeltime ! for diag
      use init_module,only: istep ! for debug
      use fields_module,only: tlbc
!
! Calculate line integrals for o2,o,n2,o3,co2,no:
!
! 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)
     |  o2,   ! molecular oxygen (mmr)
     |  o1,   ! atomic oxygen (mmr)
     |  ox,   ! atomic oxygen (o+o3) (mmr)
     |  n2,   ! molecular nitrogen (mmr)
     |  o3,   ! ozone (mmr)
     |  co2,  ! carbon dioxide (mmr)
     |  no,   ! nitric oxide (mmr)
     |  barm  ! mean molecular weight
!
! Vertical and slant-column output integrations will be column 
! number density (molecules/cm^2).
!
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) :: 
     |  vo2,  ! o2 vertical integration
     |  vo1,  ! o1 vertical integration
     |  vn2,  ! n2 vertical integration
     |  vo3,  ! o3 vertical integration
     |  vco2, ! co2 vertical integration
     |  vno,  ! no vertical integration
     |  sco2, ! o2 slant column integration
     |  sco1, ! o1 slant column integration
     |  scn2, ! n2 slant column integration
     |  sco3, ! o3 slant column integration
     |  scco2,! co2 slant column integration
     |  scno  ! no 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,dimension(nlonp4,nlat) ::       ! global (not just subdomain)
     |  sin_chi,     ! sin(chi)
     |  cos_chi,     ! cos(chi)
     |  rt_sinchi,   ! sqrt(sin(chi))
     |  slt          ! local time
      real,dimension(lev0:lev1,lon0:lon1) :: ! for diagnostic plotting
     |  chi_plt, slt_plt, idn_plt, ridn
!
! Set up sun related geometry according to current model date and time.
! Also set day/night index idn (init_module).
!
      rtpi = sqrt(pi)
      rlat = -.5*pi+(float(lat-1)+.5)*dphi
      coslat=cos(rlat)
      sinlat=sin(rlat)
      idn(lon0:lon1,lat) = 1
      do i=lon0,lon1
        slt(i,lat)=amod(secs/3600.+(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))
        if (chi(i,lat) > 1.8326) idn(i,lat) = 0

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

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

!     call addfld('CHI',' ',' ',chi_plt,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('SLT',' ',' ',slt_plt,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('IDN',' ',' ',idn_plt,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! 6/2/06 btf: addfld does not seem to handle this correctly:
!     call addfld('CHI2D',' ',' ',chi ,'lon',lon0,lon1,'lat',lat,lat,0)
!     call addfld('SLT2D',' ',' ',slt ,'lon',lon0,lon1,'lat',lat,lat,0)
!     ridn = float(idn)
!     call addfld('IDN2D',' ',' ',ridn,'lon',lon0,lon1,'lat',lat,lat,0)

!
! rp = z+re
      do k=lev0,lev1
        rp(k,:) = z(k,lon0:lon1)+re
      enddo
!     call addfld('RP',' ',' ',rp,'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! ti is tn at interfaces:
!     ti(1,:) = tn(lev1,lon0:lon1) ! tn bottom boundary is stored in top slot
      ti(1,:) = tlbc(lon0:lon1,lat) ! Lower boundary is in tlbc
      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 addfld('TNI',' ',' ',ti,'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
!     call addfld('CHAP_O2',' ',' ',o2,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('CHAP_O1',' ',' ',o1,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('CHAP_N2',' ',' ',n2,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('CHAP_O3',' ',' ',o3,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('CHAP_CO2',' ',' ',co2,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Calculate line integrals (vo[x] and sco[x] are output):
!
! o2:
      call line_integ(o2,rmass_o2,vo2,sco2,ti,barm,
     |  cos_chi(lon0:lon1,lat),sin_chi(lon0:lon1,lat),
     |  rt_sinchi(lon0:lon1,lat),rp,idn(lon0:lon1,lat),
     |  lon0,lon1,lev0,lev1,lat)
! o1:
      call line_integ(o1,rmass_o1,vo1,sco1,ti,barm,
     |  cos_chi(lon0:lon1,lat),sin_chi(lon0:lon1,lat),
     |  rt_sinchi(lon0:lon1,lat),rp,idn(lon0:lon1,lat),
     |  lon0,lon1,lev0,lev1,lat)
! n2:
      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,lat),
     |  lon0,lon1,lev0,lev1,lat)
! o3:
      call line_integ(o3,rmass_o3,vo3,sco3,ti,barm,
     |  cos_chi(lon0:lon1,lat),sin_chi(lon0:lon1,lat),
     |  rt_sinchi(lon0:lon1,lat),rp,idn(lon0:lon1,lat),
     |  lon0,lon1,lev0,lev1,lat)
! co2:
      call line_integ(co2,rmass_co2,vco2,scco2,ti,barm,
     |  cos_chi(lon0:lon1,lat),sin_chi(lon0:lon1,lat),
     |  rt_sinchi(lon0:lon1,lat),rp,idn(lon0:lon1,lat),
     |  lon0,lon1,lev0,lev1,lat)
! no:
      call line_integ(no,rmass_no,vno,scno,ti,barm,
     |  cos_chi(lon0:lon1,lat),sin_chi(lon0:lon1,lat),
     |  rt_sinchi(lon0:lon1,lat),rp,idn(lon0:lon1,lat),
     |  lon0,lon1,lev0,lev1,lat)
!
! Save column number densities:
!     call addfld('VO2',' ',' ',vo2(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('VO1',' ',' ',vo1(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('VN2',' ',' ',vn2(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('VO3',' ',' ',vo3(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('VCO2',' ',' ',vco2(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('VNO',' ',' ',vno(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Save slant line integrals:
!     call addfld('SCO2',' ',' ',sco2(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('SCO1',' ',' ',sco1(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('SCN2',' ',' ',scn2(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('SCO3',' ',' ',sco3(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('SCCO2',' ',' ',scco2(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('SCNO',' ',' ',scno(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,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)
      use params_module,only: dz
      use cons_module,only: expz,grav,p0,avo,gask,pi,re,check_exp
      use init_module,only: istep ! for debug
      implicit none
!
! Args:
      integer,intent(in) :: lon0,lon1,lev0,lev1,lat
      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:
      integer :: k,i
      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
!
!     write(6,"('Enter line_integ: lon0,1=',2i3,' lev0,1=',2i3,
!    |  ' lat=',i3,' fmin,max=',2e12.4)") lon0,lon1,lev0,lev1,lat,
!    |  minval(f),maxval(f)

! 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) then ! daytime
            if (.not.check_exp) then
              r2ig(k,i) = 2.*v(k,i)*exp(exparg)*
     |          rtpi*rt_sinchi(i)*rtrp(k,i)
            else
              r2ig(k,i) = 2.*v(k,i)*expo(exparg)*
     |          rtpi*rt_sinchi(i)*rtrp(k,i)
            endif
          else
            r2ig(k,i) = big
          endif
        enddo
      enddo
!
! Slant line integrals (0 if obscured by earth):
      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
