#include "dims.h"
      SUBROUTINE HOXCHEM
      use cons_module,only: len1,len2,len3,kmax,kmaxp1,kmaxm1,shapiro,
     |  p0,boltz,dtx2inv,expz
      implicit none
C     ****
C     ****     This subroutine deals with HOX chemistry implicitly by
C     ****      determining psistar(HOX,n+1) from:
C     ****
C     ****      psistar(HOX,n+1) = (-1./(2*Dt) + sqrt(1/(2*Dt)**2+4*
C     ****
C     ****       (Q*N*MBAR/RMTRU)(n) * (psi(HOX,n-1)*dtx2inv+P*RMTRU/
!                From cons_mod.F:  dtx2inv = 1./(2*dt)
C     ****
C     ****       (N*MBAR))(n)))/(2*Q*N*MBAR/RMTRU)(n)
C     ****
C     ****      Here S1 = P = number density source independent of
C     ****            n(HOX)
C     ****
C     ****           S2 = Q = (number density loss proportional to
C     ****            n(HOX)**2)/n(HOX)**2
C     ****
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
#include "cmpbnd.h"
#include "compcom.h"
!
! Local:
      integer :: i,k
      integer :: npjp2k,npjp1k,npnmk,npjm1k,npjm2k,nmsk,ntk,nhoxk,
     |  nhoxnmk
C     ****
C     ****     Calculate psi(HOX,n-1) in S15 using Shapiro smoother
C     ****
C     ****     Stage 1
C     ****
      NPJP2K = NJP2 + NPHOXNM
      NPJP1K = NJP1 + NPHOXNM
      NPNMK = NJ + NPHOXNM
      NPJM1K = NJM1 + NPHOXNM
      NPJM2K = NJM2 + NPHOXNM
      DO 1 I = 1,LEN3
        S14(I,1) = F(I,NPNMK)-shapiro*(F(I,NPJP2K)+F(I,NPJM2K)-4.*
     1    (F(I,NPJP1K)+F(I,NPJM1K))+6.*F(I,NPNMK))
    1 CONTINUE
C     ****
C     ****     Stage 2
C     ****
      DO 2 I = 3,LEN3-2
        S15(I,1) = S14(I,1)-shapiro*(S14(I+2,1)+S14(I-2,1)-4.*
     |    (S14(I+1,1)+S14(I-1,1))+6.*S14(I,1))
    2 CONTINUE
C     ****
C     ****     S14 = N*MBAR = p0*exp(-s)*MBAR/(k*T)
C     ****
      NMSK = NJ + NMS -1
      NTK = NJ + NT -1
      DO 3 K = 1,KMAX
        NMSK = NMSK + 1
        NTK = NTK + 1
        DO 3 I = 1,LEN1
          S14(I,K) = p0*expz(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(boltz*
     1      F(I,NTK))
    3 CONTINUE
C     ****
C     ****     S13  = psistar(HOX,n+1)
C     ****
      DO 4 I = 1,LEN2
      S12(I,1) = (dtx2inv+S3(I,1))**2+4.*S2(I,1)*S14(I,1)/
     1    RMTRU(I,1)*(S15(I,1)*dtx2inv+S1(I,1)*RMTRU(I,1)/S14(I,1))
!     S12(I,1) = merge(S12(I,1),1.E-20,S12(I,1)-1.E-20>=0.)
      if (s12(i,1) < 0.) s12(i,1) = 1.e-20
      S13(I,1) = 2.*(S15(I,1)*dtx2inv+S1(I,1)*RMTRU(I,1)/S14(I,1))/
     1    (dtx2inv+S3(I,1)+SQRT(S12(I,1)))
C
    4 CONTINUE
C     ****
C     ****     Extrapolate to level KMAXP1
C     ****
      DO 5 I = 1, LEN1
        S13(I,KMAXP1) = 2.*S13(I,KMAX)-S13(I,KMAXM1)
    5 CONTINUE
C     ****
C     ****     Place psi(HOX,n-1) in NJNPNM(HOX)
C     ****      S13 in psi(HOX,,n-1)
C     ****
      NHOXK = NJNP+NPHOXNM
      NHOXNMK = NJ+NPHOXNM
      DO 6 I = 1,LEN3
        F(I,NHOXK) = F(I,NHOXNMK)
        F(I,NHOXNMK) = S13(I,1)
    6 CONTINUE
      RETURN
      END