#include "dims.h"
      SUBROUTINE CMPHE
      use cons_module,only: len1,len3,kmax,kmaxm1,kmaxp1,imax,imaxp2,
     |  expz,grav,expzmid,gask,t0
      implicit none
C     ****     ADVANCES HE COMPOSITION BY ONE TIME STEP
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
#include "helium.h"
#include "mwt.h"
#include "diffk.h"
#include "cmpdat.h"
!
! Local:
      real, parameter :: ALFAHE=-0.38, PSHEB=0.1154E-5
      real :: PHIHE(3)=(/.270,.404,.322/)
      integer :: ibndb,ibnd,i,ntk,nmsk,nphek,nwk
      real :: xyhe
C     ****
C     ****     BOUNDARY CONDITIONS
C     ****     VALUE OF PSI GIVEN AT UPPER AND LOWER BOUNDARIES
C     ****     IBNDB = 1 (T1*DPSI/DZ + T2*PSI + T3 = 0)
C     ****     IBND = 2 (T4 = VALUE OF PSI AT UPPER BOUNDARY)
C     ****     LOWER BOUNDARY
C     ****
      IBNDB = 0
      IBND = 2
      DO 1 I=1,LEN1
      T1(I)=0.
      T2(I)=1.
      T3(I)=-SQRT(WHE(1)*WHE(2))
    1 CONTINUE
C     ****    UPPER BOUNDARY
      NTK=NJ+NT+KMAXM1
      NMSK=NJ+NMS+KMAX
      NPHEK=NJ+NPHE+KMAX
      NWK=NJ+NW+KMAX
      IF(IHE.EQ.0)THEN
C     ****     SET FOLD TO CURRENT VALUE OF MBAR*T**(3./2.)*PSH
	  DO 2 I=1,IMAX+1
	  FOLD(I,J)=F(I+2,NMSK)*(T0(KMAXP1)+F(I+2,NTK))**(3./2.)*.5*
     1    (F(I+2,NPHEK)+F(I+2,NPHEK-1))
    2     CONTINUE
C     ****     SET T4 TO CURRENT VALUE OF PSH AT UPPER BNDRY
	  DO 3 I=3,LEN1-2
	  T4(I)=.5*(F(I,NPHEK)+F(I,NPHEK-1))
    3     CONTINUE
      ELSE
	  DO 4 I=3,LEN1-2
	  T4(I)=.5*(FOLD(I-2,J)+FOLD(I-2,J+1))/(F(I,NMSK)*(T0(KMAXP1)+
     1    F(I,NTK))**(3./2.))
    4     CONTINUE
      ENDIF
C     ****     PERIODIC POINTS
      DO 5 I=1,2
      T4(I)=T4(I+IMAX)
      T4(I+IMAXP2)=T4(I+2)
    5 CONTINUE
      DO 8 I=1,IMAX+1
      PSHT(I,J)=T4(I+2)
    8 CONTINUE
C     ****     SOURCES ARE ZERO
      DO 6 I=1,LEN3
      S1(I,1)=0.
      S2(I,1)=0.
    6 CONTINUE
      XYHE = 1.E-10
      CALL MINOR(NPHE,NPHENM,RMHE,PHIHE,ALFAHE,IBND,IBNDB,WHE,XYHE,
     1  NPDHHE,difkk)
C*************************************************
C     NPHEK = NJNP + NPHE - 1
C     DO K = 1,KMAXP1,4
C       DO I = 1,IMAX
C         PLOT(I,J,K) = F(I+2,NPHEK + K)
C       ENDDO
C     ENDDO
C     IF(J.EQ.JMAX)THEN
C       DO K = 1,KMAXP1,4
C         CALL EZCNTR(PLOT(1,1,K),IMAX,JMAX)
C         WRITE(6,*)'K = ', K
C       ENDDO
C     ENDIF
C*************************************************
C     ****
C     ****     UPDATE ALFA AND BETA
C     ****     NOTE THAT:
C     ****     S5 = G
C     ****     S6 = GAMA
C     ****     S11 = AX
C     ****     S12 = EX
      DO 7 I=3,LEN1-1
      BETA(I-2,J)=9.2E-6*grav*expz(KMAX)*expzmid*(1.+T0(KMAXP1)/3330.)/
     |  (gask*S11(I,KMAXP1))
      ALFA(I-2,J)=(2./dz*(1.+S6(I,KMAX))/(1.-S6(I,KMAX))-S12(I,KMAXP1)+
     |  expz(KMAX)*expzmid*F(I,NWK)/S11(I,KMAXP1))/(F(I,NMSK)*
     |  (T0(KMAXP1)+F(I,NTK))**(3./2.)*BETA(I-2,J))
      BETA(I-2,J)=2.*S5(I,KMAX)/(dz*(1.-S6(I,KMAX))*BETA(I-2,J))
    7 CONTINUE
C     ****     PERIODIC POINTS
      ALFA(IMAX+1,J)=ALFA(1,J)
      BETA(IMAX+1,J)=BETA(1,J)
      RETURN
      END