SUBROUTINE CHAPMN
      implicit none
C     ****      CALCULATE LINE INTEGRAL FOR EACH SPECIES AT EACH GRID
C     ****      POINT IN NNO2, NNO, NNN2
      include "params.h"
      include "blnk.h"
      include "vscr.h"
      include "cons.h"
      include "index.h"
      include "buff.h"
      include "phys.h"
      include "strt.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))
      RLAT=-.5*C(110)*FLOAT(1-NHEMI)+(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

      write(6,"('chapmn: j=',i3,' secs=',f10.3,' c(4)=',f8.2,
     |  ' t1=',/(6e12.4))") j,secs,c(4),t1
      write(6,"('  t2=',/(6e12.4))") t2 
      write(6,"('  t3=',/(6e12.4))") t3
      write(6,"('  t4=',/(6e12.4))") t4 

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)
          S7(I,1)=merge(S7(I,1),0.,S7(I,1)>=0.)
    7     CONTINUE
      ELSE
          NPSK=NJ+NPS+(N-1)*KMAXP1
          DO 8 I=1,LEN3
          S7(I,1)=F(I,NPSK)
          S7(I,1)=merge(S7(I,1),0.,S7(I,1)>=0.)
    8     CONTINUE
      ENDIF
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
      S5(I,1)=S7(I,1)*C(112)*S6(I,1)*merge(0.56498823/(0.06651874+
     1S5(I,1)),(1.0606963+0.5564383*S5(I,1))/((S5(I,1)+1.7245609)*
     2S5(I,1)+1.0619896),S5(I,1)-8.>=0.)
   13 CONTINUE
C     ****      S4=2.*IG
      FACTOR=C(54)*RMASS(N)/C(57)
      DO 14 K=1,KMAXP1
      DO 14 I=1,LEN1
      S4(I,K)=2.*S7(I,K)*EXP(S9(I,K)*(1.-T2(I))*FACTOR/S8(I,K))*C(112)*
     1T4(I)*S6(I,K)
   14 CONTINUE
C     ****      SN=LINE INTEGRAL  (0.0 IF OBSCURED BY EARTH)
      NNK=NNO2+(N-1)*KMAXP1-1
      DO 15 K=1,KMAXP1
      NNK=NNK+1
      DO 15 I=1,LEN1
      F(I,NNK)=merge(1.E80,S4(I,K)-S5(I,K),T3(I)>=0.)
      F(I,NNK)=merge(F(I,NNK),1.E80,S9(I,K)*T2(I)-C(51)>=0.)
      F(I,NNK)=merge(S5(I,K),F(I,NNK),T3(I)>=0.)
   15 CONTINUE
    6 CONTINUE

      call addfsech('NNO2',f(1,nno2),zimxp,kmaxp1,kmaxp1,j)
      call addfsech('NNO' ,f(1,nno ),zimxp,kmaxp1,kmaxp1,j)
      call addfsech('NNN2',f(1,nnn2),zimxp,kmaxp1,kmaxp1,j)

      RETURN
      END
C