SUBROUTINE CHAPMN
      use init_module,only: secs
      implicit none
C     ****      CALCULATE LINE INTEGRAL FOR EACH SPECIES AT EACH GRID
C     ****      POINT IN NNO2, NNO, NNN2
      include "params.h"
      include "fcom.h"
      include "vscr.h"
      include "cons.h"
      include "index.h"
      include "buff.h"
      include "phys.h"
!
! Local:
      real :: rlat,coslat,sinlat,factor
      integer :: i,nzk,ntk,n,k,nnk,nps1k,nps2k,npsk,nmsk,kk,nnvk
C     ****      T1=SOLAR ZENITH ANGLE,CHI
C     ****      T2=SIN(CHI)
C     ****      T3=COS(CHI)
C     ****      T4=SQRT(SIN(CHI))
! 1/00: removing nhemi
!     RLAT=-.5*C(110)*FLOAT(1-NHEMI)+(FLOAT(J-1)+.5)*C(2)
      RLAT=-.5*C(110)+(FLOAT(J-1)+.5)*C(2)
      COSLAT=COS(RLAT)
      SINLAT=SIN(RLAT)
      DO 23 I=1,LEN1
      T1(I)=FLOAT(I-3)
   23 CONTINUE
      DO 1 I=1,LEN1
C     ****      T1=LOCAL TIME
      T1(I)=AMOD(SECS/3600.+(T1(I)*C(1)+C(110))*12./C(110),24.)
C     ****      T1=CHI
      T1(I)=ACOS(C(95)*SINLAT+C(96)*COSLAT*COS(C(110)*(T1(I)-12.)/12.))
      T2(I)=SIN(T1(I))
      T3(I)=COS(T1(I))
      T4(I)=SQRT(T2(I))
    1 CONTINUE
!
! idn(zimxp) is day-night index (day=1, night=0) for current latitude
!   (see fcom.h):
      idn(:) = 1
      do i=1,len1
        if (t1(i) > 1.8326) idn(i) = 0
      enddo

!     write(6,"('chapmn: j=',i2,' idn=',/,(20i3))") j,idn

C     ****      S9=RP
      NZK=NJ+NZ
      DO 2 I=1,LEN3
      S9(I,1) = F(I,NZK)+C(51)
    2 CONTINUE
C     ****      S8=T (INTERPOLATED TO KMAXP1 LEVELS)
C     ****      LEVEL 1
      NTK=NJ+NT+KMAX
      DO 3 I=1,LEN1
      S8(I,1)=F(I,NTK)+T0(1)
    3 CONTINUE
C     ****      LEVELS 2 THRU KMAX
      NTK=NJ+NT-1
      DO 4 K=2,KMAX
      NTK=NTK+1
      DO 4 I=1,LEN1
      S8(I,K)=.5*(F(I,NTK)+F(I,NTK+1))+T0(K)
    4 CONTINUE
C     ****      LEVEL KMAXP1
      NTK=NJ+NT+KMAX-1
      DO 5 I=1,LEN1
      S8(I,KMAXP1)=F(I,NTK)+T0(KMAXP1)
    5 CONTINUE
C     ****      SPECIES DO LOOP
      DO 6 N=1,3
C     ****      S7=PSI(N)
      IF(N.EQ.3)THEN
        NPS1K=NJ+NPS
        NPS2K=NPS1K+KMAXP1
        DO 7 I=1,LEN3
          S7(I,1)=1.-F(I,NPS1K)-F(I,NPS2K)
    7   CONTINUE
      ELSE
        NPSK=NJ+NPS+(N-1)*KMAXP1
        DO 8 I=1,LEN3
          S7(I,1)=F(I,NPSK)
    8   CONTINUE
      ENDIF
      do i=1,len3
        if (s7(i,1) < 0.) s7(i,1) = 0.
      enddo
C     ****      S7=NO. DENSITY INTEGRAL
C     ****      TOP OF MODEL
C     ****      S7(KMAXP1)=N0*P0*PSI(N)*EXPS*MBAR/(RMASS(N)**2*G)
      FACTOR=C(85)*C(81)*EXPS(KMAX)*C(86)/(RMASS(N)**2*C(54))
      NMSK=NJ+NMS+KMAX
      DO 9 I=1,LEN1
      S7(I,KMAXP1)=FACTOR*.5*(S7(I,KMAX)+S7(I,KMAXP1))*F(I,NMSK)
    9 CONTINUE
C     ****      INTEGRATE DOWN TO LEVEL 1
      DO 10 KK=1,KMAX
      K=KMAXP1-KK
      FACTOR=C(85)*C(81)*EXPS(K)/(RMASS(N)*C(54))*C(3)
      DO 10 I=1,LEN1
      S7(I,K)=S7(I,K+1)+FACTOR*S7(I,K)
   10 CONTINUE
C     ****      COPY VERTICAL COLUMN DENSITIES TO NNVO2, NNVO, NNVN2
      NNVK=NNVO2+(N-1)*KMAXP1
      DO 16 I=1,LEN3
      F(I,NNVK)=S7(I,1)
   16 CONTINUE
C     ****      S6=SQRT(RP/2HP)
      FACTOR=RMASS(N)*C(54)/(2.*C(57))
      DO 11 I=1,LEN3
      S6(I,1)=SQRT(S9(I,1)*FACTOR/S8(I,1))
   11 CONTINUE
C     ****      S5=YP
      DO 12 K=1,KMAXP1
      DO 12 I=1,LEN1
      S5(I,K)=S6(I,K)*ABS(T3(I))
   12 CONTINUE
C     ****      S5=IP
      DO 13 I=1,LEN3
      if (s5(i,1) >= 8.) then
        s5(i,1) = S7(I,1)*C(112)*S6(I,1)*(0.56498823/(0.06651874+
     |            S5(I,1)))
      else
        s5(i,1) = S7(I,1)*C(112)*S6(I,1)*((1.0606963+0.5564383*
     |    S5(I,1))/((S5(I,1)+1.7245609)*S5(I,1)+1.0619896))
      endif
   13 CONTINUE
C     ****      S4=2.*IG
      FACTOR=C(54)*RMASS(N)/C(57)
      DO K=1,KMAXP1
        DO I=1,LEN1
!
! FP overflow here on prospect (probably at exp(x) where x is too big)
! 3/10/00: problem at j=13. Compile with -fpe2 to avoid abort.
!
          S4(I,K)=2.*S7(I,K)*EXP(S9(I,K)*(1.-T2(I))*FACTOR/S8(I,K))*
     |      C(112)*T4(I)*S6(I,K)
        enddo
      enddo
C     ****      SN=LINE INTEGRAL  (0.0 IF OBSCURED BY EARTH)
      NNK=NNO2+(N-1)*KMAXP1-1
      DO K=1,KMAXP1
        NNK=NNK+1
        DO I=1,LEN1
          if (t3(i) >= 0.) then
            f(i,nnk) = 1.e80
          else
            f(i,nnk) = s4(i,k)-s5(i,k)
          endif
          if (S9(I,K)*T2(I)-C(51) < 0.) f(i,nnk) = 1.e80
          if (t3(i) >= 0.) f(i,nnk) = s5(i,k)
        enddo
      enddo
    6 CONTINUE
      RETURN
      END
C