#include "dims.h"
      SUBROUTINE CMPN4S
      use cons_module,only: len1,len2,kmax,kmaxp1,brn2d,expz,p0,
     |  expzmid_inv, boltz, rmassinv_o2, rmassinv_o, rmassinv_n2
      implicit none
C     ****
C     ****     ADVANCES PSI(N4S) BY ONE TIME STEP
C     ****
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "crates_const.h"
#include "crates_tdep.h"
#include "compcom.h"
#include "mwt.h"
#include "phys.h"
#include "diffk.h"
#include "cmpdat.h"
!
! Local:
      real :: PHIN4S(3)=(/0.651,0.731,0.741/)
      real :: alfa,xyn4s
      integer :: i,k,ibnd,ibndb
      integer :: npnok,npno2k,npn4sk,npn2dk,npohk,ntk,nps1k,nps2k,npo1k,
     |  nmsk,nno2k,nek,nqtefk,no2pk,nopk,nn2pk,nnpk,nnopk,nphoxk
C     ****
C     ****     BOUNDARIES
C     ****
      NPNOK = NJ+NPNO
      NPNO2K = NJ+NPNO2
      NPN4SK = NJ+NPN4S
      NPN2DK = NPN2D+NJNP
      NPOHK = NJ+NPOH
      NTK = NJ+NT+KMAX
      NPS1K = NJ+NPS
      NPS2K = NPS1K+KMAXP1
      NPO1K = NJ+NPO1
      NMSK = NJ+NMS
      NNO2K = NNO2
      NEK = NJ+NE
      NQTEFK = NQTEF
      NO2PK = NJ+NO2P
      NOPK = NJ+NOP
      NN2PK = NN2P
      NNPK = NNP
      NNOPK = NNOP
      NPHOXK = NJ+NPHOX
      DO 1 I=1,LEN1
C     ****
C     ****     T5 = N*MBAR
C     ****
      T5(I) = p0*expz(1)*expzmid_inv*F(I,NMSK)/(boltz*F(I,NTK))
C     ****     VALUE AT BOTTOM GIVEN BY PHOTOCHEMICAL EQUILIBRIUM.
      T1(I) = 0.
      T2(I) = 1.
! .5* changed to SQRT in 1st continuation line below as per modsrc.kibo:
      T3(I)= -(F(I,NQTEFK)*(1.-BRN2D)+T5(I)*(F(I,NPN2DK)/RMN2D*
     |   (T5(I)*BETA4*F(I,NPO1K)*rmassinv_o+BETA5(I,1)*SQRT(F(I,NEK)+
     |   F(I,NEK+1))+BETA7)+.5*(XJNO(I,1)+XJNO(I+LEN1,1))*F(I,NPNOK)/
     |   RMNO)+T5(I)*(RK2(I,1)*F(I,NOPK)*(1.-F(I,NPS1K)-F(I,NPO1K))*
     |   rmassinv_n2+RK6*F(I,NNPK)*F(I,NPS1K)*rmassinv_o2+RK8*
     |   F(I,NNPK)*F(I,NPO1K)*rmassinv_o)+SQRT(F(I,NEK)*F(I,NEK+1))
     |   *(RA1(I,1)*F(I,NNOPK)*0.15+RA3(I,1)*F(I,NN2PK)*1.1))/
     |   (T5(I)*(BETA1(I,1)*F(I,NPS1K)*rmassinv_o2+BETA3(I,1)*
     |   F(I,NPNOK)/RMNO+BETA8*F(I,NPHOXK)/RMTRU(I,1)*RATIO1(I,1)*
     |   RATIO3(I,1)+BETA13*F(I,NPNO2K)/RMNO2)-RK4*F(I,NO2PK))*
     |   RMN4S/T5(I)
C     ****
C     ****     ZERO DIFFUSIVE FLUX AT TOP
C     ****
      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     ****
C     ****     S3 = EXPS,  S2 = RN4S
C     ****
      S2(I,K) = .5*(F(I,NQTEFK)+F(I,NQTEFK+1))*(1.-BRN2D)
    2 CONTINUE
      NTK = NJ+NT
      DO 3 I=1,LEN2
C     ****
C     ****     S3 = N*MBAR  (K+1/2)
C     ****
      S3(I,1)= XNMBAR(I,1)
C     ****
C     ****     S2 = NUMBER DENSITY PRODUCTION INDEPENDENT OF N(N4S)
C     ****
! .5* changed to SQRT in 1st continuation line below as per modsrc.kibo:
      S2(I,1)= S2(I,1)+S3(I,1)*(F(I,NPN2DK)/RMN2D*(S3(I,1)*BETA4*
     |  F(I,NPO1K)*rmassinv_o+BETA5(I,1)*SQRT(F(I,NEK)+F(I,NEK+1))+
     |  BETA7)+.5*(XJNO(I,1)+XJNO(I+LEN1,1))*F(I,NPNOK)/RMNO)+S3(I,1)
     |  *(RK2(I,1)*F(I,NOPK)*(1.-F(I,NPS1K)-F(I,NPO1K))*rmassinv_n2
     |  +RK6*F(I,NNPK)*F(I,NPS1K)*rmassinv_o2+RK8*F(I,NNPK)*F(I,NPO1K)*
     |  rmassinv_o)+SQRT(F(I,NEK)*F(I,NEK+1))*(RA1(I,1)*F(I,NNOPK)*0.15+
     |  RA3(I,1)*F(I,NN2PK)*1.1)
C     ****
C     ****     S1 = (NUMBER DENSITY LOSS)/N(N4S)
C     ****
      S1(I,1) = -S3(I,1)*(BETA1(I,1)*F(I,NPS1K)*rmassinv_o2+BETA3(I,1)*
     |  F(I,NPNOK)/RMNO+BETA8*F(I,NPHOXK)/RMTRU(I,1)*RATIO1(I,1)*
     |  RATIO3(I,1)+BETA13*F(I,NPNO2K)/RMNO2)-RK4*F(I,NO2PK)
    3 CONTINUE
      IBND = 0
      IBNDB = 0
      ALFA = 0.
      XYN4S = 1.E-6
      CALL MINOR(NPN4S,NN4SNM,RMN4S,PHIN4S,ALFA,IBND,IBNDB,WN4S,
     1XYN4S,NPDHN4S,difkk)
      RETURN
      END
C