#include "dims.h"
!
      subroutine chapmn(idn)
      use init_module,only: secs,iday
      implicit none
!
! Time-gcm wrapper for chapman.
!
#include "params.h"
#include "fcom.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
!
! Args:
      integer,intent(out) :: idn(zimxp)
!
! Local:
      logical,save :: isfirst=.true.
!
      if (isfirst) then
!       write(6,"('CHAPMN: Calling tgcm/waccm chapman routine.')")
        isfirst = .false.
      endif
      call chapman(secs,iday,
     |  idn ,f(1,nj+nps),f(1,nj+npo1),f(1,nj+npo3),f(1,nj+nps2),
     |  f(1,nj+nt),f(1,nj+nz) ,f(1,nj+nms), zimxp,zkmxp,zjmx,j,
     |  f(1,nnvo2),f(1,nnvo),f(1,nnvn2),f(1,nnvo3),f(1,nno2),f(1,nno),
     |  f(1,nnn2),f(1,nno3))
!
! Update f-array day/night index from chapman's idn.
! (note there is only one level in f(i,ndn)).
!
      f(:,ndn) = float(idn(:))
      end subroutine chapmn
!
!-----------------------------------------------------------------------
!
      subroutine chapman(secs,iday,idn,xo2mmr,xommr,xo3mmr,xoxmmr,
     |                   t,z,rmbar,nlon,nlev,nlat,lat,vo2,vo,vn2,vo3,
     |                   sco2,sco,scn2,sco3)
!
!-----------------------------------------------------------------------
! Calculate vertical and slant-column integrals for O2, O, N2 on the
! geographic grid.
!
! Calling Sequence:
!   tgcm: Call from wrapper chapmn using f-array (fcom.h), see above.
!         Chapmn is called from dynamics.F.
!
! Code history:
!   Standardized for tiegcm and waccm by Ben Foster, August 2000.
!
! Notes:
!   - N2 mass mixing ratio is calculated from input as 1.-xo2mmr-xommr.
!   - F90 array syntax is used only in the longitude dimension.
!   - tgcm and waccm vertical indices are reversed (tgcm k=1 == ccm k=plev).
!   - Input nlon includes periodic points.
!   - Input nlev includes extra level (zkmxp, plev).
!   - Day/night index rdn(nlon) is set here, and will be used in qrj.
!   - This version is for time-gcm only. Chapman for tiegcm does not have
!     ozone.
!-----------------------------------------------------------------------
!
! Parameters and constants from init module:
      use cons_module,only: expz,dphi,dlamda,avo,p0,re,grav,dz,pi,
     |                      rmass_o2,rmass_o,rmass_n2,rmass_o3,gask
      implicit none
!
! Input arguments:
      integer,intent(in) :: 
     |  iday,                ! current calendar day
     |  lat,                 ! current latitude index
     |  nlon,                ! number of longitudes      (zimxp, plon)
     |  nlev,                ! number of pressure levels (zkmxp, plev)
     |  nlat                 ! number of latitudes       (zjmx , plat)
      real,intent(in) :: 
     |  secs,                ! current seconds in the day
     |  xo2mmr(nlon,nlev),   ! O2 mass mixing ratio
     |  xommr (nlon,nlev),   ! O  mass mixing ratio
     |  xo3mmr(nlon,nlev),   ! O3 mass mixing ratio
     |  xoxmmr(nlon,nlev),   ! OX mass mixing ratio (O1+O3)
     |  t     (nlon,nlev),   ! neutral temperature (deg K)
     |  z     (nlon,nlev),   ! geopotential height (cm)
     |  rmbar (nlon,nlev)    ! mbar mean mass (1/(o2/32+o/16+n2/28))
!
! Output arguments:
      integer,intent(out) :: idn(nlon) ! day/night index at current lat
      real,intent(out) :: 
     |  vo2 (nlon,nlev),     ! vertical O2 column density (cm3)
     |  vo  (nlon,nlev),     ! vertical O  column density (cm3)
     |  vn2 (nlon,nlev),     ! vertical N2 column density (cm3)
     |  vo3 (nlon,nlev),     ! vertical O3 column density (cm3)
     |  sco2(nlon,nlev),     ! slant column (line integral) O2 (cm3)
     |  sco (nlon,nlev),     ! slant column (line integral) O  (cm3)
     |  scn2(nlon,nlev),     ! slant column (line integral) N2 (cm3)
     |  sco3(nlon,nlev)      ! slant column (line integral) O3 (cm3)
!
! Local:
      real,parameter :: big=1.e80
      real :: rlat,coslat,sinlat,factor,x,exparg,delta,rtpi
      real :: fmin,fmax,fmin0,fmax0
      real :: 
     |  chi(nlon),            ! solar zenith angle
     |  sinchi(nlon),         ! sin(chi)
     |  coschi(nlon),         ! cos(chi)
     |  rtsinchi(nlon)        ! sqrt(sin(chi))
      real ::
     |  rp(nlon,nlev),        ! z+re         (S9)
     |  tn(nlon,nlev),        ! t interpolated to interfaces (S8)
     |  rtrp(nlon,nlev),      ! sqrt(rp/2hp) (S6)
     |  yp(nlon,nlev),        ! yp and ip    (S5)
     |  r2ig(nlon,nlev)       ! 2.*ig        (S4)
      real :: 
     |  vcomp(nlon,nlev,4),   ! vertical integrals for o2,o,n2,o3
     |  sccomp(nlon,nlev,4),  ! slant-line integrals for o2,o,n2,o3
     |  rmass(4)              ! molecular weights o2,o,n2,o3
      integer :: i,k,kk,n
      real :: slt(nlon)
      integer :: old_fpe_flags, new_fpe_flags
!
! Exec:
!
      delta = atan(tan(23.5*pi/180.)*sin(2.*pi*float(iday-80)/365.))
      rlat = -.5*pi+(float(lat-1)+.5)*dphi
      coslat=cos(rlat)
      sinlat=sin(rlat)
      rtpi = sqrt(pi)  ! C(112)
!
! Set up sun related geometry according to current model date and time:
      do i=1,nlon
        chi(i)=float(i-3)            ! T1
        chi(i)=amod(secs/3600.+(chi(i)*dlamda+pi)*12./pi,24.) ! local time
        slt(i) = chi(i)
        chi(i)=acos(sin(delta)*sinlat+cos(delta)*coslat*
     |    cos(pi*(chi(i)-12.)/12.))
        sinchi(i)=sin(chi(i))        ! T2
        coschi(i)=cos(chi(i))        ! T3
        rtsinchi(i)=sqrt(sinchi(i))  ! T4
      enddo
!
! Day/night index rdn(nlon) (this will be used in qrj):
      idn(:) = 1 ! init all day
      do i=1,nlon
        if (chi(i) >= 1.8326) idn(i) = 0 ! night
      enddo
!
! rp (S9):
      do k=1,nlev
        rp(:,k) = z(:,k)+re
      enddo
!
! Input temperature is at midpoints. Interpolate tn (S8) to interfaces:
! (no extrapolation)
      tn(:,1) = t(:,nlev) ! tgcm bottom boundary tn is stored in top slot
      do k=2,nlev-1
        tn(:,k) = .5*(t(:,k-1)+t(:,k)) 
      enddo
      tn(:,nlev) = t(:,nlev-1)
!
! Loop over o2,o,n2,o3:
      do n=1,4
!
! Initialize vertical columns with input mass mixing ratios:
        select case(n)
          case(1)
            rmass(n) = rmass_o2
            do k=1,nlev
              vcomp(:,k,n) = xo2mmr(:,k)
            enddo
          case(2)
            rmass(n) = rmass_o
            do k=1,nlev
              vcomp(:,k,n) = xommr (:,k)
            enddo
          case(3)
            rmass(n) = rmass_n2
            do k=1,nlev
!
! Use N2 = 1.-O2-OX, where OX = O1+O3.
              vcomp(:,k,n) = 1.-xo2mmr(:,k)-xoxmmr(:,k)
            enddo
          case(4)
            rmass(n) = rmass_o3
            do k=1,nlev
              vcomp(:,k,n) = xo3mmr(:,k)
            enddo
          case default
            write(6,"('>>> bad n=',i5)") n
            stop 'chapman'
        endselect
        where(vcomp(:,:,n) < 0.) vcomp(:,:,n) = 0. ! insure positive
!
! Top (tgcm k==nlev, ccm k==1) in cm3:
        factor = avo*p0*expz(nlev-1)*exp(-.5*dz)/(rmass(n)**2*grav)
        vcomp(:,nlev,n) = factor*.5*(vcomp(:,nlev-1,n)+vcomp(:,nlev,n))*
     |                    rmbar(:,nlev)
!
! Integrate from top down:
        do k=1,nlev-1
          kk = nlev-k ! e.g., nlev==29, kk==28,1
          factor = avo*p0*expz(kk)/(rmass(n)*grav)*dz
          vcomp(:,kk,n) = vcomp(:,kk+1,n)+factor*vcomp(:,kk,n) ! S7
        enddo
!
! Transfer vertical integrals to output arrays:
        select case(n)
          case (1)
            vo2(:,:) = vcomp(:,:,n)
!           call fminmax(vo2,nlon*nlat,fmin,fmax)
!           write(6,"('chapmn: lat=',i2,' vo2 min,max=',2e12.4)")
!    |        lat,fmin,fmax
          case (2)  
            vo (:,:) = vcomp(:,:,n)
!           call fminmax(vo,nlon*nlat,fmin,fmax)
!           write(6,"('chapmn: lat=',i2,' vo min,max=',2e12.4)")
!    |        lat,fmin,fmax
          case (3)
            vn2(:,:) = vcomp(:,:,n)
!           call fminmax(vn2,nlon*nlat,fmin,fmax)
!           write(6,"('chapmn: lat=',i2,' vn2 min,max=',2e12.4)")
!    |        lat,fmin,fmax
          case (4)
            vo3(:,:) = vcomp(:,:,n)
!           call fminmax(vo3,nlon*nlat,fmin,fmax)
!           write(6,"('chapmn: lat=',i2,' vo3 min,max=',2e12.4)")
!    |        lat,fmin,fmax
        end select
!
! Set up for slant-line integral:
        factor = rmass(n)*grav/(2.*gask)
        do k=1,nlev
          rtrp(:,k) = sqrt(rp(:,k)*factor/tn(:,k)) ! S6
          yp(:,k) = rtrp(:,k)*abs(coschi(:))       ! S5
        enddo
        do k=1,nlev
          do i=1,nlon
            if (yp(i,k) >= 8.) then
              yp(i,k) = vcomp(i,k,n)*rtpi*rtrp(i,k)*
     |          (0.56498823/(0.06651874+yp(i,k)))
            else
              yp(i,k) = vcomp(i,k,n)*rtpi*rtrp(i,k)*
     |          ((1.0606963+0.5564383*yp(i,k))/((yp(i,k)+1.7245609)*
     |            yp(i,k)+1.0619896))
            endif
          enddo ! i=1,nlon
        enddo ! k=1,nlev
        factor = grav*rmass(n)/gask
!
! The exp() in the r2ig calculation causes overflow fpe at night that 
! can be flushed to huge. However, overflows here will cause large 
! number of underflows in qrj.
!
! On compaq prospect, f90 -fpe2 allows continuation with overflows
! and underflows, but is slow due to fpe trapping. Using the below 
! conditional on idn (do calculation only in daytime) and exparg 
! getting too large (near the terminator) prevents overflows on compaq 
! and consequent underflows in qrj.
!
        do k=1,nlev
          do i=1,nlon
            exparg = rp(i,k)*(1.-sinchi(i))*factor/tn(i,k)
            if (idn(i)==1.and.exparg < 650.) then ! daytime
              r2ig(i,k) = 2.*vcomp(i,k,n)*exp(exparg)*
     |          rtpi*rtsinchi(i)*rtrp(i,k)
            else
              r2ig(i,k) = big
            endif
          enddo
        enddo
! 
! Slant line integrals (0 if obscured by earth):
        do k=1,nlev
          do i=1,nlon
            if (coschi(i) >= 0.) then
              sccomp(i,k,n) = big
            else
              sccomp(i,k,n) = r2ig(i,k)-yp(i,k)
            endif
            if (rp(i,k)*sinchi(i)-re < 0.) sccomp(i,k,n) = big
            if (coschi(i) >= 0.) sccomp(i,k,n) = yp(i,k)
          enddo
        enddo
      enddo ! n=1,4 for o2,o,n2,o3
!
! Transfer line integrals to output arrays:
!     subroutine addfsech(fname,f,idim1,idim2,ndim2,ilat)
!
      sco2(:,:) = sccomp(:,:,1) 
!     call fminmax(sco2,nlon*nlev,fmin,fmax)
!     write(6,"('chapmn: lat=',i2,' sco2 min,max=',2e12.4)")
!    |  lat,fmin,fmax
!     call addfsech('SCO2',' ',' ',sco2,nlon,nlev,nlev,lat)
!
      sco (:,:) = sccomp(:,:,2) 
!     call fminmax(sco,nlon*nlev,fmin,fmax)
!     write(6,"('chapmn: lat=',i2,' sco min,max=',2e12.4)")
!    |  lat,fmin,fmax
!     call addfsech('SCO',' ',' ',sco,nlon,nlev,nlev,lat)
!
      scn2(:,:) = sccomp(:,:,3) 
!     call fminmax(scn2,nlon*nlev,fmin,fmax)
!     write(6,"('chapmn: lat=',i2,' scn2 min,max=',2e12.4)")
!    |  lat,fmin,fmax
!     call addfsech('SCN2',' ',' ',scn2,nlon,nlev,nlev,lat)
!
      sco3(:,:) = sccomp(:,:,4) 
!     call fminmax(sco3,nlon*nlev,fmin,fmax)
!     write(6,"('chapmn: lat=',i2,' sco3 min,max=',2e12.4)")
!    |  lat,fmin,fmax
!     call addfsech('SCO3',' ',' ',sco3,nlon,nlev,nlev,lat)
!
!     new_fpe_flags = for_set_fpe(old_fpe_flags)
      return
      end subroutine chapman