#include "dims.h"
      SUBROUTINE CMPN4S
      use cons_module,only: len1,len2,brn2d,expz,rmassinv,kmax,kmaxp1,
     |  p0,expzmid_inv,boltz
      use crates_module,only: rk8,rk4,rk6,beta4,beta7
      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_tdep.h"
#include "mwt.h"
#include "phys.h"
!
! 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
      real :: fn4s(zimxp,zkmxp)
!
      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)=p0*expz(1)*expzmid_inv*F(I,NMSK)/(boltz*F(I,NTK))

!       write(6,"('comp_n4s: lat=',i2,' i=',i2,' barm=',e12.4,
!    |    ' xnmbar=',e12.4,' tn=',e12.4)") 
!    |    j,i,f(i,nmsk),t5(i),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*T5(I)*F(I,NPS2K)*rmassinv(2)+BETA5(I,1)*
     2        F(I,NEK)+BETA7)+BETA8(I,1)*F(I,NPNOK)/RMNO)
     3        /(BETA1(I,1)*F(I,NPS1K)*rmassinv(1)+BETA3(I,1)*
     4        F(I,NPNOK)/RMNO+T5(I)*BETA17(I,1)*F(I,NPS2K)*
     5        rmassinv(2)*(1.-F(I,NPS2K)-F(I,NPS1K))*rmassinv(3))

!       write(6,"('comp_n4s: lat=',i2,' i=',i2)") j,i
!       write(6,"('  xnmbar_local(lev0,i)=',e12.4)") t5(i)
!       write(6,"('  qtef (lev0,i) =',e12.4)") f(i,nqtefk)
!       write(6,"('  n2d  (lev0,i) =',e12.4)") f(i,npn2dk)
!       write(6,"('  o1   (lev0,i) =',e12.4)") f(i,nps2k)
!       write(6,"('  beta5(lev0,i) =',e12.4)") beta5(i,1)

!       write(6,"('  n4s_out(lev0,i)=',e12.4)") t3(i)

C     ****     ZERO DIFFUSIVE FLUX AT TOP
      T4(I)=0.
    1 CONTINUE

!     write(6,"('comp_n4s: lat=',i2,' n4s lbc=',/,(4e15.7))")
!    |  j,t3(:)

!     do k=1,kmax
!       fn4s(:,k) = t3(:)
!     enddo
!     call addfsech('N4S_OUT',' ',' ',fn4s,zimxp,zkmxp,zkmx,j)

C     ****
C     ****     SOURCES
C     ****
      NQTEFK=NQTEF-1
      DO 2 K=1,KMAX
      NQTEFK=NQTEFK+1
      DO 2 I=1,LEN1
C     ****     S3 = expz,  S2 = RN4S
      S3(I,K)=expz(K)
      S2(I,K)=.5*(F(I,NQTEFK)+F(I,NQTEFK+1))*(1.-BRN2D)
    2 CONTINUE

!     call addfsech('N4S_PROD',' ',' ',s2,zimxp,zkmxp,zkmx,j)

      NTK=NJ+NT
      DO 3 I=1,LEN2

C     ****     S3 = N*MBAR  (K+1/2)
      S3(I,1)=S3(I,1)*p0/(boltz*F(I,NTK)*(F(I,NPS1K)*rmassinv(1)+
     1F(I,NPS2K)*rmassinv(2)+(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3)))

C     ****     S2 = NUMBER DENSITY PRODUCTION INDEPENDENT OF N(N4S)
      S2(I,1)=S2(I,1)+S3(I,1)*(F(I,NPN2DK)/RMN2D*(S3(I,1)*BETA4*
     |  F(I,NPS2K)*rmassinv(2)+BETA5(I,1)*.5*(F(I,NEK)+F(I,NEK+1))+
     |  BETA7)+.5*(BETA8(I,1)+BETA8(I+LEN1,1))*F(I,NPNOK)/RMNO)+
     |  S3(I,1)*(RK2(I,1)*F(I,NOPK)*(1.-F(I,NPS1K)-F(I,NPS2K))*
     |  rmassinv(3)+RK6*F(I,NNPK)*F(I,NPS1K)*rmassinv(1)+RK8*
     |  F(I,NNPK)*F(I,NPS2K)*rmassinv(2))+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     ****     S1 = (NUMBER DENSITY LOSS)/N(N4S)
      S1(I,1)=-S3(I,1)*(BETA1(I,1)*F(I,NPS1K)*rmassinv(1)+BETA3(I,1)*
     1F(I,NPNOK)/RMNO+S3(I,1)*BETA17(I,1)*F(I,NPS2K)*rmassinv(2)*(1.
     2-F(I,NPS2K)-F(I,NPS1K))*rmassinv(3))-RK4*F(I,NO2PK)
    3 CONTINUE

!     call addfsech('XNMBAR',' ',' ',s3,zimxp,zkmxp,zkmx,j)
!     call addfsech('N4S_PROD',' ',' ',s2,zimxp,zkmxp,zkmx,j)
!     call addfsech('N4S_LOSS',' ',' ',s1,zimxp,zkmxp,zkmx,j)

      IBND=0
      IBNDB=0
      ALFA=0.
      CALL MINOR(NPN4S,NN4SNM,RMN4S,PHIN4S,ALFA,IBND,IBNDB,0)

!     call addfsech('N4S_OUT',' ',' ',f(1,njnp+npn4s),zimxp,zkmxp,
!    |  zkmx,j)
!     call addfsech('N4S_TM1',' ',' ',f(1,njnp+nn4snm),zimxp,zkmxp,
!    |  zkmx,j)

      RETURN
      END
C