SUBROUTINE CMPN4S
      implicit none
C     ****
C     ****     ADVANCES PSI(N4S) BY ONE TIME STEP
C     ****
      include "params.h"
      include "blnk.h"
      include "vscr.h"
      include "cons.h"
      include "index.h"
      include "buff.h"
      include "crates.h"
      real :: rmn4s,rmn2d,rmno,brn2d,cee
      COMMON/MASS/ RMN4S,RMN2D,RMNO,BRN2D,CEE
!
! Local:
      real :: PHIN4S(3),alfa
      integer :: npnok,npn4sk,npn2dk,ntk,nps1k,nps2k,nmsk,nno2k,nek,
     |  nqtefk,no2pk,nopk,nn2pk,nnpk,nnopk,i,k,ibnd,ibndb
!
      DATA RMN4S,RMN2D,RMNO/14.,14.,30./
      DATA PHIN4S/0.651,0.731,0.741/
C     ****
C     ****     BOUNDARIES
C     ****
      NPNOK=NJ+NPNO
      NPN4SK=NJ+NPN4S
      NPN2DK=NPN2D+NJNP
      NTK=NJ+NT+KMAX
      NPS1K=NJ+NPS
      NPS2K=NPS1K+KMAXP1
      NMSK=NJ+NMS
      NNO2K=NNO2
      NEK=NJ+NE
      NQTEFK=NQTEF
      NO2PK=NJ+NO2P
      NOPK=NJ+NOP
      NN2PK=NN2P
      NNPK=NNP
      NNOPK=NNOP
      DO 1 I=1,LEN1
C     ****     T5 = N*MBAR
      T5(I)=C(81)*EXPS(1)*C(87)*F(I,NMSK)/(C(84)*F(I,NTK))
C     ****     VALUE AT BOTTOM GIVEN BY PHOTOCHEMICAL EQUILIBRIUM.
      T1(I)=0.
      T2(I)=1.
      T3(I) = -RMN4S/T5(I)*(F(I,NQTEFK)*(1.-BRN2D)/T5(I)+
     1       F(I,NPN2DK)/RMN2D*(
     1        BETA4(I,1)*T5(I)*F(I,NPS2K)/RMASS(2)+BETA5(I,1)*
     2        F(I,NEK)+BETA7(I,1))+BETA8(I,1)*F(I,NPNOK)/RMNO)
     3        /(BETA1(I,1)*F(I,NPS1K)/RMASS(1)+BETA3(I,1)*
     4        F(I,NPNOK)/RMNO+T5(I)*BETA17(I,1)*F(I,NPS2K)/
     5        RMASS(2)*(1.-F(I,NPS2K)-F(I,NPS1K))/RMASS(3))
C     ****     ZERO DIFFUSIVE FLUX AT TOP
      T4(I)=0.
    1 CONTINUE
C     ****
C     ****     SOURCES
C     ****
      NQTEFK=NQTEF-1
      DO 2 K=1,KMAX
      NQTEFK=NQTEFK+1
      DO 2 I=1,LEN1
C     ****     S3 = EXPS,  S2 = RN4S
      S3(I,K)=EXPS(K)
      S2(I,K)=.5*(F(I,NQTEFK)+F(I,NQTEFK+1))*(1.-BRN2D)
    2 CONTINUE
      NTK=NJ+NT
      DO 3 I=1,LEN2
C     ****     S3 = N*MBAR  (K+1/2)
      S3(I,1)=S3(I,1)*C(81)/(C(84)*F(I,NTK)*(F(I,NPS1K)/RMASS(1)+
     1F(I,NPS2K)/RMASS(2)+(1.-F(I,NPS1K)-F(I,NPS2K))/RMASS(3)))
C     ****     S2 = NUMBER DENSITY PRODUCTION INDEPENDENT OF N(N4S)
      S2(I,1)=S2(I,1)+S3(I,1)*(F(I,NPN2DK)/RMN2D*(
     1S3(I,1)*BETA4(I,1)*F(I,NPS2K)/RMASS(2)+BETA5(I,1)*.5*(F(I,NEK)+
     2F(I,NEK+1))+BETA7(I,1))+.5*(BETA8(I,1)+BETA8(I+LEN1,1))*F(I,NPNOK)
     3/RMNO)+S3(I,1)*(RK2(I,1)*F(I,NOPK)*(1.-F(I,NPS1K)-F(I,NPS2K))/
     4RMASS(3)+RK6(I,1)*F(I,NNPK)*F(I,NPS1K)/RMASS(1)+RK8(I,1)*
     5F(I,NNPK)*F(I,NPS2K)/RMASS(2))+SQRT(F(I,NEK)*F(I,NEK+1))*(RA1(I,1)
     6*F(I,NNOPK)*0.15+RA3(I,1)*F(I,NN2PK)*1.1)
C     ****     S1 = (NUMBER DENSITY LOSS)/N(N4S)
      S1(I,1)=-S3(I,1)*(BETA1(I,1)*F(I,NPS1K)/RMASS(1)+BETA3(I,1)*
     1F(I,NPNOK)/RMNO+S3(I,1)*BETA17(I,1)*F(I,NPS2K)/RMASS(2)*(1.
     2-F(I,NPS2K)-F(I,NPS1K))/RMASS(3))-RK4(I,1)*F(I,NO2PK)
    3 CONTINUE
      IBND=0
      IBNDB=0
      ALFA=0.
      CALL MINOR(NPN4S,NN4SNM,RMN4S,PHIN4S,ALFA,IBND,IBNDB)
      RETURN
      END
C