C
C PURPOSE: BUILD STATISTICAL EQUILIBRIUM MATRIX
C
C INPUTS: NDIM (pretty sure this is an output!! SG)
C
C OUTPUTS:
C
C COMMON:
C
C COMMENTS: OCTOBER 6, 1999, P. JUDGE
C 
	SUBROUTINE SE0_BUILD(NDIM)

C  THIS SUBROUTINES BUILDS THE COEFFICIENT MATRIX OF THE STATISTICAL
C  EQUILIBRIUM EQUATIONS IN THE B-FRAME, ONCE THE B-DIRECTION IN THE
C  S-FRAME, AND THE HEIGHT OVER THE LIMB OF THE SCATTERER ARE GIVEN.
C  A CYLINDRICALLY SYMMETRIC, UNPOLARIZED RADIATION FIELD IS ASSUMED.
C  QUANTUM COHERENCES ARE NEGLECTED.

	INCLUDE 'PREC'
	INCLUDE 'PARAM'
	INCLUDE 'CATOM'
	INCLUDE 'CSE'
	INCLUDE 'CINTS'
C added next 6jul22
	INCLUDE 'CINPUT'
	INCLUDE 'CLU'
C       
C	DIMENSION DUM(MJTOT,MJTOT)
	LOGICAL LEQK
C       
C       CALCULATE TENSORS OF RADIATION FIELD FOR ALL POSSIBLE TRANSITIONS
C       
C	WRITE (LOUT,*) 'RADIATION TENSORS'
C       
C	CALL CPTIME('SE0_BUILD',0,0,1)
C	CALL CPTIME('SE0:FIELDB',0,0,1)
	DO IJ=1,NK
	   DO IJ1=1,NK
C       
C       ENUMERATE RADIATIVELY
C       PERMITTED TRANSITIONS
C       
	      KR=KRAD(IJ,IJ1)				
	      IF (KR.NE.0) THEN 
		 CALL FIELDB(KR)
	      END IF
	   END DO
	END DO
C	CALL CPTIME('SE0:FIELDB',0,1,1)
C       
C       BUILD COEFFICIENT MATRIX  SCOEFF(I,I1)
C       
C	CALL CPTIME('SE0:R&C',0,0,1)
	I=0
	IZ=0
	DO IJ=1,NK
	   CALL R1(IJ,R1COEFF)
	   C1TMP0=CE(IJ,IZ)
	   CALL C2(IJ,C2COEFF)
	   CALL C6(NK,IJ,C6COEFF)
C Note stepping by 2 added with 6jul22 version
	   DO K=0,NINT(TWO*AJ(IJ)),2
C       
C       ROW NUMBER
C       
	      I=I+1
	      IF (K.NE.0) THEN
		 WEIGHT(I)=ZERO
	      ELSE
		 WEIGHT(I)=SQRT(TWO*AJ(IJ)+ONE)
	      END IF			
	      C1TMPK=CE(IJ,K)
	      C1COEFF=C1TMP0-C1TMPK
	      I1=0
	      DO IJ1=1,NK
		 KR=KRAD(IJ,IJ1)
		 CALL R4(IJ,IJ1,K,R4COEFF)
		 IF (IJ1.LT.IJ) CALL C3(IJ,IJ1,K,C3COEFF)
		 IF (IJ1.GT.IJ) CALL C5(IJ,IJ1,K,C5COEFF)
C Note stepping by 2 added with 6jul22 version
		 DO K1=0,NINT(TWO*AJ(IJ1)),2
		    LEQK=K1.EQ.K
C       
C       COLUMN NUMBER
C       
		    I1=I1+1					
		    SCOEFF(I,I1)=ZERO
C       
C       R(1,2,6), C(1,2,6),
C       -> IJ1=IJ
C       
		    IF (IJ1.EQ.IJ) THEN			
		       IZ=0
		       CALL R2(IJ,K,K1,IZ,R2COEFF)
		       CALL R6(NK,IJ,K,K1,IZ,R6COEFF)
		       SCOEFF(I,I1)=-R2COEFF-R6COEFF+SCOEFF(I,I1)
		       IF (LEQK) SCOEFF(I,I1)=-R1COEFF-C1COEFF-C2COEFF-C6COEFF
     * 		    +SCOEFF(I,I1)
      
		    END IF
		    IF (IJ1.LT.IJ) THEN
		       IF (KR.NE.0) THEN
			  IZ=0
			  CALL R3(KR,IJ,IJ1,K,K1,IZ,R3COEFF)
			  SCOEFF(I,I1)=R3COEFF+SCOEFF(I,I1)
		       END IF
		       IF (LEQK) SCOEFF(I,I1)=C3COEFF+SCOEFF(I,I1)
		    END IF
		    IF (IJ1.GT.IJ) THEN
		       IF (KR.NE.0) THEN
			  IZ=0
			  CALL R5(KR,IJ,IJ1,K,K1,IZ,R5COEFF)
			  SCOEFF(I,I1)=R5COEFF+SCOEFF(I,I1)
		       END IF
		       IF (LEQK) SCOEFF(I,I1)=R4COEFF+C5COEFF+SCOEFF(I,I1)
		    END IF
		 END DO
	      END DO
	   END DO
	END DO
C debug statement added 6jul22
	if(IDEBUG .GT. 0) WRITE(LOUT,*)' SE0_BUILD MADE A  ',I,' X ',I1, 'MATRIX'
	IF (I.NE.I1) CALL STOP('SE_BUILD: SCOEFF IS NOT SQUARE')
	NDIM=I
C	CALL CPTIME('SE0:R&C',0,1,1)
C	CALL CPTIME('SE0_BUILD',0,1,1)
	RETURN
	END