#include "dims.h"
      SUBROUTINE CMPN2D
      use cons_module,only: len1,len2,kmax,kmaxp1,brn2d,rmassinv,
     |  expz,p0,boltz
      use crates_module,only: beta7,rk10,beta6,beta2,beta4
      implicit none
C     ****
C     ****     CALCULATES PSI(N2D) ASSUMING PHOTOCHEMICAL EQUILIBRIUM
C     ****     NOTE: THIS IS DIAGNOSTIC.  CALL BEFORE CMPN4S AND CMPNO
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:
      integer :: npnok,npn4sk,npn2dk,ntk,nps1k,nps2k,nno2k,nek,
     |  nqtefk,nopk,nn2pk,nnopk,k,i
      real :: xnmbar(zimxp,zkmxp)
!
!     call addfsech('QTEF'  ,' ',' ',f(1,nqtef),zimxp,zkmxp,zkmx,j)
!     call addfsech('N2P'   ,' ',' ',f(1,nn2p ),zimxp,zkmxp,zkmx,j)
!     call addfsech('NOP'   ,' ',' ',f(1,nnop ),zimxp,zkmxp,zkmx,j)
!     call addfsech('O1_N2D',' ',' ',f(1,nj+nps2),zimxp,zkmxp,zkmx,j)
!     call addfsech('NE_N2D',' ',' ',f(1,nj+ne)  ,zimxp,zkmxp,zkmx,j)
!     call addfsech('RA1'   ,' ',' ',ra1         ,zimxp,zkmxp,zkmx,j)
!     call addfsech('RA3'   ,' ',' ',ra3         ,zimxp,zkmxp,zkmx,j)
!     call addfsech('RK3'   ,' ',' ',rk3         ,zimxp,zkmxp,zkmx,j)

      NPNOK=NJ+NPNO      ! no
      NPN4SK=NJ+NPN4S    ! n4s
      NPN2DK=NPN2D+NJNP  ! n2d (itc)
      NTK=NJ+NT          ! tn
      NPS1K=NJ+NPS       ! o2
      NPS2K=NPS1K+KMAXP1 ! o1
      NNO2K=NNO2         ! no2
      NEK=NJ+NE          ! ne
      NQTEFK=NQTEF-1     ! qtef (qrj)
      NOPK=NJ+NOP        ! op
      NN2PK=NN2P         ! n2p
      NNOPK=NNOP         ! nop
      DO 1 K=1,KMAX
      NQTEFK=NQTEFK+1
      DO 1 I=1,LEN1
C     ****     S1 = expz(K+1/2),    S2 = PN2D(K+1/2)
      S1(I,K)=expz(K)
      S2(I,K)=.5*(F(I,NQTEFK)+F(I,NQTEFK+1))*BRN2D
    1 CONTINUE
      DO 3 I=1,LEN2
C     ****     S1 = N*MBAR  (K+1/2)
      S1(I,1)=S1(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)))
      xnmbar(i,1) = s1(i,1)
C     ****     S2 = TOTAL N2D PRODUCTION
      S2(I,1)=S2(I,1)+RK3(I,1)*F(I,NN2PK)*S1(I,1)*F(I,NPS2K)*
     |  rmassinv(2)+(RA1(I,1)*F(I,NNOPK)*0.85+RA3(I,1)*F(I,NN2PK)*0.9)*
     |  SQRT(F(I,NEK)*F(I,NEK+1))
C     ****     S3 = (TOTAL N2D LOSS)/N(N2D)
      S3(I,1)=S1(I,1)*(BETA2*F(I,NPS1K)*rmassinv(1)+BETA4*
     1F(I,NPS2K)*rmassinv(2)+BETA6*F(I,NPNOK)/RMNO)+BETA7+
     2BETA5(I,1)*SQRT(F(I,NEK)*F(I,NEK+1))+RK10*F(I,NOPK)
C     ****     S1 = PSI(N2D)
      S1(I,1)=RMN2D*S2(I,1)/(S3(I,1)*S1(I,1))
    3 CONTINUE

!     call addfsech('XNMBARM' ,' ',' ',xnmbar,zimxp,zkmxp,zkmx,j)
!     call addfsech('N2D_PROD',' ',' ',s2    ,zimxp,zkmxp,zkmx,j)
!     call addfsech('N2D_LOSS',' ',' ',s3    ,zimxp,zkmxp,zkmx,j)

      DO 4 I=1,LEN2
C     ****     F(NPN2DK) = PSI(N2D)
      F(I,NPN2DK)=S1(I,1)
    4 CONTINUE

!     call addfsech('N2D_UPD',' ',' ',f(1,npn2d+njnp),
!    |  zimxp,zkmxp,zkmx,j)

      RETURN
      END
C