c c cGENMIX c
Encyclopaedia Index c

c************************************************************************
c   A message to prospective users of the GENXIX program
c
c   From: Brian Spalding  (dbs@cham.co.uk)
c   Date: 28.07.99
c   Text: This program is being made available for down-loading from the
c         CHAM web-site in response to requests from prospective users.
c
c         If there is sufficient interest, I can make relevant parts of 
c         the text of the original book available also. Please let me
c         know if you would like this.
c
c************************************************************************
c                    gnmx.for                06.10.97
c
C   This Fortran file was recreated during 1997 from printed records
c   of the 1977 coding. It appears to work correctly.
c
c   How to use it is explained in the book entitled:
c
c  'GENMIX, a general computer program for two-dimensional
c   parabolic phenomena', by D. B. Spalding               
c
c   published by Pergamon Press, London, 1977
c
c************************************************************************
c
      BLOCK DATA
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
C     'GENMIX, A GENERAL COMPUTER PROGRAM FOR TWO-DIMENSIONAL
C      PARABOLIC PHENOMENA', BY D. B. SPALDING
C     REPORT NO. HTS/77/9, FEBRUARY 1977,
C     IMPERIAL COLLEGE, MECHANICAL ENGINEERING DEPARTMENT, LONDON,SW72BX
C
C      APPENDIX A (BASIC PROGRAM) - COMBUSTION OF METHANE AND AIR IN
C                  A DIVERGENT DUCT EXHAUSTING INTO THE ATMOSPHERE.
C-----------------------------------------------------------------------
CHAPTER 1  1  1  1  1  1  1  1  PRELIMINARIES  1  1  1  1  1  1  1  1
      COMMON/COMA/
     1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
     2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
     3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
     4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
     5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
     6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
      COMMON/COMB/
     1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
     2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
     3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
     4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
     5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
     6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
     7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
     8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
     9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
     1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
     2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
      common/dbg/debug
      logical debug
C   
C      ------------------------------------------------------------------
      save
      DATA KUDIF/-1/
      DATA KASE,IRUN/0,0/
      DATA ITEST/1/
      DATA BIG,TINY/1.E20,1.E-10/
ccc      data big,tiny/1.e5,1.e-5/
C   
C-----------------------------------------------------------------------
CHAPTER 2  2  2  2  2  2  2  GRID AND GEOMETRY  2  2  2  2  2  2  2  2
C     ------------------------------ GRID
      DATA N,OMPOW/20,1./
c      ----- set idimf= dimension for i, ensure that idimf.ge.n
      DATA IDIMF/20/
c      ----- set krad=1 for plane, 2 for axial, 3 for point symmetry
      DATA KRAD,CSALFA/2,1./
C      ------------------------------ GEOMETRY
      DATA HEX0,XHEX0,AHEX,BHEX,CHEX/.05,0.,.01,0.,0./
      DATA HIN0,XHIN0,AHIN,BHIN,CHIN/.02,0.,0.,0.,0./
      DATA HDIV/.025/
      DATA XU,XEND,XOUT/0.,.5,1./
      DATA LASTEP,XULAST/100,2./
C      ----- FOR KIND=3, XU, XHEX0 AND XHIN0=.25 (MAIN CH.2)
C   
C-----------------------------------------------------------------------
chapter 3  3  3  3  3  3  3  dependent variables  3  3  3  3  3  3  3
c      ----- set nf= number of dependent variables, excluding velocity
       DATA NF/3/,JH,JP,JF,JOX,JTE,JPR/1,2,3,4,5,6/
c     -----set novel=1 for no velocity, novel=2 otherwise
       DATA NOVEL/2/
C
C-----------------------------------------------------------------------
CHAPTER 4  4  4  4  4  4  4  4  PROPERTY DATA  4  4  4  4  4  4  4  4
C     -------------------- S.I. UNITS
       DATA AGRAV,GASCON/9.81,8314./
c     ----- set model=1 for laminar flow,
c     ----- set model=2 for 'mixing-length' model of turbulence
       DATA MODEL/2/
ccc       data model/1/
       DATA AK,FR,CEBU,EWALL/.435,.033,.4,9./
c     ----- initial almg, and adjusted value for ch.7
       DATA ALMG/.09/, ALMGD/.075,.1,.14,.075/
c     --- set inert=1 for inert fluid, inert=2 for chemically reactive
       DATA INERT/2/
c     ------------------------------ materials
c     ---------- thermodynamics
       DATA CFU,COX,CPR,CMIX/4*1100./
ccc       data cfu,cox,cpr,cmix/4*1./
       DATA WFU,WOX,WPR,WMIX/16.,32.,28.,29./
       DATA HFU/4.E7/
ccc       data hfu/0.0/
C     ---------- CHEMICAL
       DATA STOICH,ARRCON,PREEXP/4.,18.E3,1./
C     ---------- TRANSPORT
       DATA VISFU,VISOX,VISPR,VISMIX/4*1.E-6/
       DATA PRLAM,PRTURB/.7,.86/
       DATA H,UFAC/.9,.01/
C
C-----------------------------------------------------------------------
CHAPTER 5  5  5  5  5  5  5  STARTING VALUES  5  5  5  5  5  5  5  5  
       DATA PRESS/1.E5/
C     ----- STREAM B IS PURE FUEL
       DATA UB,TB,FUB,OXB/100.,1000.,1.,0./
ccc       data ub,tb,fub,oxb/100.,100.,1.,0./
C     ----- STREAM C IS AIR
       DATA UC,TC,FUC,OXC/50.,1000.,0.,.232/
ccc       data uc,tc,fuc,oxc/100.,100.,0.,.232/
c     ----- set kex and kin for initial boundary type,
c     ----- 1 for wall, 2 for free boundary, 3 for symmetry axis
       DATA KEX,KIN/1,1/
C
C-----------------------------------------------------------------------
CHAPTER 6  6  6  6  6  6  6  6 STEP CONTROL 6  6  6  6  6  6  6  6  6
       DATA FRA,DXMAX,DXRAT/1.,1.,5./DXPSI,DXRE,DXY/3*0.0/
C     -----ENTRAINMENT CONTROL
       DATA ULIM,PEILIM,FACEXP/.02,.05,.5/
C     -----STARTING VALUES
       DATA FACE,FACI,RATE,RATI/4*1./
C
C-----------------------------------------------------------------------
CHAPTER 7  7  7  7  7  7  7  7  BOUNDARY CONDITIONS  7  7  7  7  7  7
C     ---------- STREAM A, THROUGH CENTRAL PIPE
        DATA UA,TA,FUA,OXA/100.,2000.,0.,0./
ccc        data ua,ta,fua,oxa/100.,100.,0.,0./
C     ---------- STREAM D, SURROUNDING ATMOSPHERE
        DATA TD,FUD,OXD/300.,0.,.232/
ccc        data td,fud,oxd/100.,0.,.232/
C     ----- UD IS SUPPLIED BY WAY OF THE UEX FUNCTION
C     ---------- VELOCITY ALONG OUTER BOUNDARY
        DATA UEX0,XUEX0,AUEX,BUEX,CUEX/5*0./
C     ---------- WALL TEMPERATURE OF OUTER TUBE
        DATA TWALL/299./
ccc        data twall/100./
C
C-----------------------------------------------------------------------
CHAPTER 11 11 11 11 11 11 11 11 11 11  PRINT  11  11  11  11  11  11
c     --- set ilplot=2 for  down-stream plot, =1 for no plot
c     --- set itplot=2 for cross-stream plot, =1 for no plot
        DATA ILPLOT,ITPLOT/2,2/
c     --- set nstat, nprof, nplot to number of steps between output of
c     --- station values, profiles and cross-stream plots respectively
        DATA NSTAT,NPROF,NPLOT/12,12,10000/
ccc        data nstat,nprof,nplot/1,1,10000/
c     ----- after xu=xout, nstat and nprof are set =24 at main, ch.11
C
      END
c-----------------------------------------------------------------------      
      PROGRAM MAIN
      common/dbg/debug
      logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
C-----------------------------------------------------------------------
CHAPTER 1  1  1  1  1  1  1  1  PRELIMINARIES  1  1  1  1  1  1  1  1
C-----------------------------------------------------------------------
      COMMON/COMA/
     1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
     2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
     3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
     4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
     5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
     6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
       COMMON/COMB/
     1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
     2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
     3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
     4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
     5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
     6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
     7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
     8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
     9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
     1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
     2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
      save
C
C     ------------------------------------------------------------------
C     ----- FUNCTIONS FOR BOUNDARY CONDITIONS
      HEX(X)=HEX0+X*(AHEX+X*(BHEX+X*CHEX))
      HIN(X)=HIN0+X*(AHIN+X*(BHIN+X*CHIN))
      UEX(X)=UEX0+X*(AUEX+X*(BUEX+X*CUEX))
      debug=.false.
C
C-----------------------------------------------------------------------
CHAPTER 2  2  2  2  2  2  2  GRID AND GEOMETRY  2  2  2  2  2  2  2  2 
C     SEE DATA
c     ----- kind is an index which denotes a particular geometry type
      if(debug) write(6,*)'>>> main'
      KIND=4
      IF(KRAD.EQ.1) KIND=2
      IF(KRAD.EQ.2 .AND. nint(CSALFA).EQ.1) KIND=1
      IF(KRAD.EQ.2 .AND. nint(CSALFA).EQ.0) KIND=3
C       ----- MODIFICATIONS TO DATA
      IF(KIND.NE.3) GO TO 21
      XU=.25
      XHEX0=.25
      XHIN0=.25
   21	CONTINUE
C
      SNALFA=SQRT(1.-CSALFA**2)
C       ----- STARTING VALUES
c!!!! change needed to avoid values too large for conversion
cccc  IEND=IFIX(XEND*1.E6)
cccc  IOUT=IFIX(XOUT*1.E6)
      IEND=IFIX(XEND*1.E4)
      IOUT=IFIX(XOUT*1.E4)
C       -------------------- SUBROUTINE COMP UTE, ENTRY INIT
      CALL INIT
C       ------------------------------ GRID
      DO 20 I=1,N
   20	OM(I)=(FLOAT(I-1)/FLOAT(NM1))**OMPOW
C
C     -------------------- SUBROUTINE COMP UTE, ENTRY GRID
      CALL GRID
C
C-----------------------------------------------------------------------
CHAPTER 3  3  3  3  3  3  3  DEPENDENT VARIABLES  3  3  3  3  3  3  3
C     SEE DATA
C     U(I)= VELOCITY
C     F(I,JH)= STAGNATION ENTHALPY
C     F(I,JP)= PHI= OXIDANT CONCENTRATION - F(I,JF)*STOICH
C     F(I,JF)= FUEL CONCENTRATION
C     F(I,JOX)= OXIDANT CONCENTRATION
C     F(I,JTE)= TEMPERATURE
C     F(I,JPR)= PRODUCT CONCENTRATION
C
C-----------------------------------------------------------------------
CHAPTER 4  4  4  4  4  4  4  4  PROPERTY DATA  4  4  4  4  4  4  4  4
C     SEE DATA
      RECWFU=1./WFU
      RECWOX=1./WOX
      RECWPR=1./WPR
      RECWMX=1./WMIX
      GAMMA=CMIX/(CMIX-GASCON*RECWMX)
      gamma=1.4
      DO 40 J=1,NF
        PRL(J)=PRLAM
        RECPRL(J)=1./PRLAM
   40 RECPRT(J)=1./PRTURB
C
C-----------------------------------------------------------------------
CHAPTER 5  5  5  5  5  5  5  STARTING VALUES  5  5  5  5  5  5  5  5
C     SEE DATA
      WB=1./(FUB*RECWFU+OXB*RECWOX+(1.-FUB-OXB)*RECWPR)
      RHOB=PRESS*WB/(TB*GASCON)
      WC=1./(FUC*RECWFU+OXC*RECWOX+(1.-FUC-OXC)*RECWPR)
      RHOC=PRESS*WC/(TC*GASCON)
      FLOB=RHOB*UB*(HDIV-HIN0)
      FLOC=RHOC*UC*(HEX0-HDIV)
      IF(KRAD.EQ.1) GO TO 55
      XSIN=XU*SNALFA
      HCOS=.5*CSALFA
      FLOB=FLOB*(XSIN+HCOS*(HDIV+HIN0))
      FLOC=FLOC*(XSIN+HCOS*(HEX0+HDIV))
   55	CONTINUE
      OMDIV=FLOB/(FLOB+FLOC+TINY)
      TMIN=.5*AMIN1(TA,TB,TC,TD,TWALL)
      tmin=90.
C     -------------------- SEQUENCE TO PUT CELL BOUNDARY AT OMDIV.
      IF(OMDIV.LE.1.E-10.OR.OMDIV.GE.(1.-1.E-10)) GO TO 53
      DO 52 I=3,NM1
      IF(OMINT(I)-OMDIV) 52,53,57
   57	IDIV=I+1
      GO TO 58
   52	CONTINUE
   58	FAC=OMDIV/OMINT(IDIV-1)
      DO 59 I=2,IDIV
   59	OM(I)=OM(I)*FAC
C     -------------------- SURBROUTINE COMP UTE, ENTRY GRID
      CALL GRID
   53	CONTINUE
C     ---------- INSERTION INTO ARRAYS
      ENTHB=TB*(CFU*FUB+COX*OXB+CPR*(1.-FUB-OXB))+.5*UB**2+HFU*FUB
      ENTHC=TC*(CFU*FUC+COX*OXC+CPR*(1.-FUC-OXC))+.5*UC**2+HFU*FUC
      enthb=tb*(cfu*fub+cox*oxb+cpr*(1.-fub-oxb))+hfu*fub
      enthc=tc*(cfu*fuc+cox*oxc+cpr*(1.-fuc-oxc))+hfu*fuc
      PHIB=OXB-FUB*STOICH
      PHIC=OXC-FUC*STOICH
      DO 501 I=1,N
      IF(OM(I).GT.OMDIV) GO TO 503
      U(I)=UB
      F(I,JH)=ENTHB
      F(I,JP)=PHIB
      F(I,JF)=FUB
      GO TO 501
  503	U(I)=UC
      F(I,JH)=ENTHC
      F(I,JP)=PHIC
      F(I,JF)=FUC
  501	CONTINUE
C
c     --------------------------------- enter main loop at chapter 7
      GO TO 700
C
C-----------------------------------------------------------------------
CHAPTER 6  6  6  6  6  6  6  6  STEP CONTROL  6  6  6  6  6  6  6  6
C     SEE DATA
  600 DXY=FRA*Y(NM2)
      DXRE=DXY*PEI/(.5*(R(1)+R(N))*EMU(1)+TINY)
      DXINC=DXLAST*DXRAT
C
C     -------------------- DETERMINATION OF BOUNDARY TYPE
C     -------------------- I BOUNDARY
      IF(ISTEP.GE.IEND) GO TO 610
      KIN=1
      GO TO 611
  610	IF(PSII.LE.TINY) GO TO 612
      KIN=2
      GO TO 611
  612	KIN=3
C     -------------------- E BOUNDARY
  611 IF(ISTEP.GE.IOUT) GO TO 613
      KEX=1
      GO TO 614
  613 KEX=2
  614 CONTINUE   
C
C     -------------------- ENTRAINMENT RATES
      IF(KIN.NE.2.AND.KEX.NE.2) GO TO 602
      KUDIF=ISTEP
      UMAX=U(1)
      UMIN=U(1)
      DO 615 I=2,N
      UMAX=AMAX1(UMAX,U(I))
  615	UMIN=AMIN1(UMIN,U(I))
      UDIF=UMAX-UMIN
C     --------------------------- I BOUNDARY
      IF(KIN.NE.2) GO TO 601
      RATI=ABS((U(2)-U(1))/(UDIF*ULIM+TINY))
      RMI=(R(2)+R(3))*(EMU(2)+EMU(3))*RECYDF(2)*RATI/(1.+RATI)
      FACI=FACI*RATI**FACEXP
      FACI=AMAX1(0.1,AMIN1(FACI,10.))
      RMI=RMI*FACI
      IF(MODEL.EQ.2) RMI=AMAX1(RMI,0.4*UDIF*RHO(1)*R(1))
C     --------------------------- E BOUNDARY
  601	IF(KEX.NE.2) GO TO 602
      RATE=ABS((U(NM1)-U(N))/(UDIF*ULIM+TINY))
      RME=-(R(NM2)+R(NM1))*(EMU(NM2)+EMU(NM1))*RECYDF(NM2)*RATE/
     1                                                    (1.+RATE)
      FACE=FACE*RATE**FACEXP
      FACE=AMAX1(.01,AMIN1(FACE,10.))
      RME=RME*FACE
      IF(MODEL.EQ.2) RME=AMAX1(RME,-0.4*UDIF*RHO(N)*R(N))
C
  602	DXPSI=PEI*PEILIM/(RMI-RME+TINY)
C
C     --------------- SET VALUE OF DX
      DX=AMIN1(DXY,DXRE,DXINC,DXPSI,DXMAX)
C
      IF(ISTEP.GE.IEND) GO TO 605
      IF(DX.LT.(XEND-XU)) GO TO 605
c     ----- reset dx so that xu will exactly equal xend at next step
      DX=XEND-XU
      IEND=ISTEP+1
      JUSTIN=ISTEP+1
C
  605	IF(ISTEP.GE.IOUT) GO TO 606
      IF(DX.LT.(XOUT-XU)) GO TO 606
c     ----- reset dx so that xu will exactly equal xout at next step
      DX=XOUT-XU
      IOUT=ISTEP+1
      JUSTEX=ISTEP+1
C
  606	IF(PSII.GT.RMI*DX) GO TO 607
      IF(PSII.LE.TINY) GO TO 607
c     ----- reset dx so that axis is reached at next step
      DX=PSII/RMI
      JUSTIN=ISTEP+1
C
C     ----- RESET DX SO THAT XU WILL NOT EXCEED XULAST
  607	DX=AMIN1(DX,XULAST-XU)
C
C     ---------- TRAP ZERO OR NEGATIVE DX
      IF(DX.GT.0.) GO TO 608
      IFIN=2
	GO TO 1100
C
C     -----DETERMINE XD
  608	XD=XU+DX
      DXLAST=DX
C
C     ----- IF CSALFA VARIES
C           RECALCULATE IT, AND SNALFA AND HCOS, HERE, FOR X=XD
      GO TO 70
C
C-----------------------------------------------------------------------
CHAPTER 7  7  7  7  7  7  7  7  BOUNDARY CONDITIONS  7  7  7  7  7  7
  700 ASSIGN 751 TO ISTART
C     ---------- GENERAL BOUNDARY CONDITION INFORMATION
C     ---------- STREAM A,THROUGH CENTRAL PIPE
C     SEE DATA
      ENTHA=TA*(CFU*FUA+COX*OXA+CPR*(1.-FUA-OXA))+.5*UA**2+HFU*FUA
      PHIA=OXA-FUA*STOICH
      RHOA=PRESS*WPR/(TA*GASCON+TINY)
      FLOA=RHOA*UA*HIN(XEND)
      IF(KRAD.EQ.2) FLOA=FLOA*(XEND*SNALFA+HCOS*HIN(XEND))
      PSII=FLOA
      PEI=FLOB+FLOC
      PSIE=PSII+PEI
C     ---------- STREAM D, SURROUNDING ATMOSPHERE
C     SEE DATA
      XUEX0=XOUT
      UD=UEX0
      XD=XU
      ENTHD=TD*(CFU*FUD+COX*OXD+CPR*(1.-FUD-OXD))+.5*UD**2+HFU*FUD
      PHID=OXD-FUD*STOICH
C     ---------- OTHER RELATED INFORMATION
      HDUCID=HIN0
      ADUCTD=HEX0-HDUCID
      IF(KRAD.EQ.2) ADUCTD=ADUCTD*(XSIN+HCOS*(HEX0+HDUCID))
      AFLOWD=ADUCTD
   70	CONTINUE
C     --------------- BOUNDARY CONDITION FOR FORWARD STEP
C     ------------------------------------------- I BOUNDARY 
      IF(KIN-2) 731,732,733
C     ------------------------------  WALL
  731	IF(ISTEP.GT.JUSTIN) GO TO 734
      U(1)=0.
      TAUI=0.
      RMI=0.
      DO 735 J=1,NF
      IBIN(J)=2
  735 RJTOTI(J)=0.    
C     -------------------- ADJUST INNER HEIGHT
  734 HIND=HIN(XD-XHIN0)
      GO TO 740
C     ------------------------------ FREE BOUNDARY
  732 IF(ISTEP.GT.JUSTIN) GO TO 736
      TAUI=0.
      U(1)=UA
      RHO(1)=RHOA
      RECRU(1)=1./(RHO(1)*U(1)+TINY)
      F(1,JH)=ENTHA
      F(1,JP)=PHIA
      F(1,JF)=FUA
      AREA=HDUCID
      IF(KRAD.EQ.2) AREA=AREA*(XU*SNALFA+HCOS*HDUCID)
      AFLOWD=AFLOWD+AREA
c!!!! mistyping
cccc      IF(ISTEP.GT.0) GO TO 740
      IF(ISTEP.eq.0) GO TO 740
  736	U(1)=U(1)+DX*AGRAV*(RHO(N)-RHO(1))*RECRU(1)
      GO TO 740
C     ------------------------------ SYMMETRY AXIS
  733	IF(ISTEP.GT.JUSTIN) GO TO 740
      TAUI=0.
      RMI=0.
      PSII=0.
      HIND=0.
      U(1)=U(2)
      DO 737 J=1,NF
  737	F(1,J)=F(2,J)
C     ---------- NO SUBSEQUENT CHANGE NEEDED
  740	CONTINUE
C
C     ------------------------------------------- E BOUNDARY
      IF(KEX-2) 741,742,743
C     ------------------------------ WALL
  741	IF(ISTEP.GT.JUSTEX) GO TO 744
C     ---------- FIRST STEP ONLY
      U(N)=0.
      RME=0.
      TAUE=0.
      IBEX(JH)=1
      F(N,JOX)=OXC
      F(N,JPR)=1.-OXC-FUC
      DO 745 J=2,NF
        IBEX(J)=2
  745 RJTOTE(J)=0.
C     ----- ADJUST ENTHANPY FIT COMPOSITION
  744	CMIX=CFU*F(N,JF)+COX*F(N,JOX)+CPR*F(N,JPR)
      F(N,JTE)=TWALL
      F(N,JH)=CMIX*F(N,JTE)+F(N,JF)*HFU
C     -------------------- ADJUST EXTERNAL HEIGHT
      HEXD=HEX(XD-XHEX0)
      GO TO 750
C     -------------------------------------------- FREE BOUNDARY
  742	IF(ISTEP.GT.JUSTEX) GO TO 746
      F(N,JH)=ENTHD
      F(N,JP)=PHID
      F(N,JF)=FUD
      F(N,JOX)=OXD
      F(N,JPR)=1.-F(N,JF)-F(N,JOX)
      VMIX=F(N,JF)*RECWFU+F(N,JOX)*RECWOX+F(N,JPR)*RECWPR
      F(N,JTE)=TD
      RHO(N)=PRESS/(VMIX*F(N,JTE)*GASCON)
      U(N)=UD
      RECRU(N)=1./(RHO(N)*U(N)+TINY)
C     ---------- ADJUSTMENT OF MIXING LENGTH CONSTANT
      ALMG=ALMGD(KIND)
C     ---------- ADJUSTMENT OF DOWNSTREAM VELOCITY
  746	UD=UEX(XD-XUEX0)
      GO TO 750
C     ---------- NO SYMMETRY AXIS
  743	CONTINUE
  750	GO TO ISTART, (751,800)
  751	ASSIGN 800 TO ISTART
      GO TO 900
C
C-----------------------------------------------------------------------
CHAPTER 8  8  8  8  8  8  8  8  ADVANCE  8  8  8  8  8  8  8  8  8  8
C     ------------------------  MOMENTUM SOURCES
C     -------------------------------------- PRESSURE GRADIENT
  800 continue
      IF(KEX.NE.2) GO TO 821
      DP=(U(N)-UD)/RECRU(N)
      GO TO 823
  821	AFLOWU=AFLOWD
      HDUCID=0.
      IF(KIN.EQ.1) HDUCID=HIND
      ADUCTD=HEXD-HDUCID
      IF(KRAD.EQ.2) ADUCTD=ADUCTD*(XD*SNALFA+HCOS*(HEXD+HDUCID))
      DA=ADUCTD-AFLOWU
      DP=DA/DADP
C     -------------------- WALL SHEAR AND MASS ADDITION
      UBAR=0.
      DO 824 I=2,NM1
  824	UBAR=UBAR+(BOM(I)*U(I))
      IF(KIN.EQ.2) UBAR=(UBAR-U(1))*PEI/PSIE+U(1)
      UBAR=(UBAR-U(1))*PEI/PSIE+U(1)
      DP=DP+DX*(-TAUI*R(1)-TAUE*R(N)+2.*RME*UBAR)/ADUCTD
      DP=AMIN1(DP,.5*DPMAX)
C
C     ------------------------------ COMP
  823 CALL SOLVE
C
C-----------------------------------------------------------------------
CHAPTER 9  9  9  9  9  9  9  9  COMPLETE  9  9  9  9  9  9  9  9  9  9
  900 CONTINUE
C     ---------------------- IGNITION SEQUENCE
      IF(ISTEP.GT.5) GO TO 931
      IF(INERT.EQ.1) GO TO 931
      T2=.5/STOICH
      DO 932 I=2,NM1
      F(I,JF)=T2*(ABS(F(I,JP))-F(I,JP))
  932	F(I,JOX)=F(I,JP)+STOICH*F(I,JF)
  931	CONTINUE
C
C     ------------------------------ THERMODYNAMIC PROPERTIES
      PRESS=PRESS+DP
      PDGSCN=PRESS/GASCON
      DO 907 I=1,N
      F(I,JOX)=AMAX1(0.,F(I,JP)+STOICH*F(I,JF))
      F(I,JPR)=1.-F(I,JF)-F(I,JOX)
      ENTH=F(I,JH)-.5*U(I)**2-HFU*F(I,JF)
      ENTH=F(I,JH)-HFU*F(I,JF)
      F(I,JTE)=ENTH/(CFU*F(I,JF)+COX*F(I,JOX)+CPR*F(I,JPR))
      IF(F(I,JTE).GT.TMIN) GO TO 941
      IF(I.EQ.1.OR.I.EQ.N) GO TO 941
      WRITE(6,942) F(I,JTE),I,ISTEP,TMIN
  942 FORMAT(27H *** TEMPERATURE, F(I,JTE)=,1PE10.3,6H AT I=,I4,7H ISTEP
     1=,I5/17H *** RESET =TMIN=,E10.3,23H *** MAIN CH.9 COMPLETE)
      F(I,JTE)=TMIN
c      write(6,*)'CFU,F(I,JF)',
c     1           CFU,F(I,JF)
c      write(6,*)'COX,F(I,JOX)',
c     1           COX,F(I,JOX)
c      write(6,*)'CPR,F(I,JPR)',
c     1           CPR,F(I,JPR)
c      write(6,*)'hfu,enth',
c     1           hfu,enth
c      write(6,*)'F(I,JH),U(I)',F(I,JH),U(I)
      stop
  941	VMIX=F(I,JF)*RECWFU+F(I,JOX)*RECWOX+F(I,JPR)*RECWPR
  907	RHO(I)=PDGSCN/(F(I,JTE)*VMIX)
      IF(KEX.EQ.1) F(N,JTE)=TWALL
      DPDX=DP/DX
C
C     ---------------------------- RADII AND Y@S
      IF(KRAD-2) 901,902,903
C     ----- KRAD=1, PLANE
  901	IF(KIN.EQ.2) HIND=ABS(PSII*RECRU(1))
      GO TO 909
C     ----- KRAD=2, AXIAL
  902	IF(KIN.NE.2) GO TO 908
      HIND=ABS(PSII*RECRU(1))
      HIND=2.*HIND/
     1         (XD*SNALFA+SQRT((XD*SNALFA)**2+2.*HIND*CSALFA)+TINY)
      GO TO 908
C     ----- KRAD=3, POINT SYMMETRY
  903	R(I)=0.
C     ----- CHANGE ABOVE STATEMENT IF NECESSARY FOR KRAD=3
      GO TO 909
  908	R(1)=XU*SNALFA+HIND*CSALFA
C     ------------------------------- COMP
  909	CALL DISTAN
C
C-----------------------------------------------------------------------
CHAPTER 10  10  10  10  10  10  10  ADJUST  10  10  10  10  10  10  10  
C
      IF(KEX.EQ.2) GO TO 1022
      AFLOWD=Y(N)+HIND-HDUCID
      IF(KRAD.EQ.2) AFLOWD=AFLOWD*(XU*SNALFA+HCOS*(Y(N)+HIND+HDUCID))
      DA1=ADUCTD/AFLOWD-1.
C       -------------------- DEPENDENCE OF AREA ON PRESSURE
      RECGMP=1./(GAMMA*PRESS)
      DADP=0.
      IF(KIN.EQ.2) DADP=PSII*RECRU(1)*(RECRU(1)*RECRU(1)*RHO(1)-RECGMP)
      SUM=0.
      DPMAX=BIG
      DO 1025 I=2,NM1
      DPMAX=AMIN1(DPMAX,RHO(I)*U(I)**2)
 1025	SUM=SUM+BOM(I)*RECRU(I)*(RECRU(I)*RECRU(I)*RHO(I)-RECGMP)
      DADP=DADP+PEI*SUM
c     -------------------- adjustment of ps,us, etc.
      IF(ABS(DA1).LT.1.E-3) GO TO 1022
      DP=DA1*AFLOWD/DADP
      DP=AMIN1(DP,.0*DPMAX)
      PRESS=PRESS+DP
      DPDX=DPDX+DP/DX
      RHOFAC=1.+DP*RECGMP
      DO 1027 I=2,NM1
      U(I)=U(I)-DP*RECRU(I)
 1027	RHO(I)=RHO(I)*RHOFAC
      IF(KIN.NE.2) GO TO 1029
      U(1)=U(1)-DP*RECRU(1)
 1029	RHO(1)=RHO(1)*RHOFAC
      RHO(N)=RHO(N)*RHOFAC
      RECRU(1)=1./(RHO(1)*U(1)+TINY)
      IF(KIN.NE.2) GO TO 1026
      HIND=ABS(PSII*RECRU(1))
      IF(KRAD.EQ.1) GO TO 1026
      HIND=2.*HIND/
     1       (XD*SNALFA+SQRT((XD*CSALFA)**2+2.*HIND*CSALFA)+TINY)
      R(I)=XU*SNALFA+HIND*CSALFA
 1026 CALL DISTAN
      AFLOWD=Y(N)+HIND-HDUCID
      IF(KRAD.EQ.2) AFLOWD=AFLOWD*(XU*SNALFA+HCOS*(Y(N)+HIND+HDUCID))
      DA2=ADUCTD/AFLOWD-1.
 1022	CONTINUE
C-----------------------------------------------------------------------
CHAPTER 11 11 11 11 11 11 11 11 11 11  PRINT  11  11  11  11  11  11 
C     SEE DATA
 1100	CONTINUE
      IF(XU.LE.XOUT) GO TO 1101
      NSTAT=24
      NPROF=24
 1101	CONTINUE     
      CALL OUTPUT
C
C-----------------------------------------------------------------------
CHAPTER 12  12  12  12  12  12  12  DECIDE  12  12  12  12  12  12  12
      IF(ISTEP.EQ.LASTEP) GO TO 1203
      IF(XU.LT.XULAST) GO TO 1202
 1203	IFIN=2
      CALL OUTPUT
 1202	IF(IFIN.EQ.1) GO TO 600
      STOP
      END

      SUBROUTINE WALL(I1,OUT1,OUT2)
C/FEB. 1977 --------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
      common/dbg/debug
      logical debug
      DIMENSION S1(2),S2(2),S3(2),S4(2),S5(2),S6(2)
      COMMON/COMA/
     1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
     2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
     3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
     4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
     5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
     6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
       COMMON/COMB/
     1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
     2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
     3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
     4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
     5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
     6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
     7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
     8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
     9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
     1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
     2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
C     
C     EFFECTS OF PRESSURE GRADIENT AND MASS TRANSFER ARE INCLUDED
C     EFFECTS OF RADIUS VARIATION ARE NEGLECTED
C     FOR VELOCITY,  OUT1=BP,     OUT2=T
C     FOR F 'S,      OUT1=FIDIF,  OUT2=T
C
CHAPTER A ------------------------------- PRELIMINARIES ----------------
      save
      DATA SHALF/.04/, BPLAST/1./
      if(debug) write(6,*)'>>> wall'
      KWALL=2-1/I1
      I2=I1+3-2*KWALL
      I3=I1+6-4*KWALL
      IF(J.GT.0) GO TO 200
CHAPTER B ------------------------------- VELOCITY ---------------------
      UREF=U(I2)
      RHOREF=RHO(I2)
      RUREF=RHOREF*UREF
      RREF=R(I2)
      VREF=EMU(I1)
      YREF=YI+(YE-YI)*OM(I1)
      RE=RUREF*YREF/VREF
      RRUREF=RREF*RUREF
      AM=(RMI-(RME+RMI)*OM(I1))/RRUREF
      EF=YREF*DP/(DX*RUREF*UREF)
      IF(MODEL.EQ.1) GO TO 110
      IF(RE.LT.132.25) GO TO 110
C     ------------------------------ TURBULENT FLOW
C     ---------- EXTENDED LOG LAW
      ER=RE*EWALL
      ARGMIN=11.5*EWALL
      NIT=0
  101 SHALF1=SHALF
      S=SHALF**2
      SLOC=S+AM+EF
      IF(SLOC.GT.0.) GO TO 104
      SLOC=TINY
      SHALF=SQRT(ABS(AM+EF))
  104 BEE=SQRT(SLOC)/AK
      ARG=ER*(SHALF+(AM/(1.+BEE)+.5*EF)/SHALF)
      IF(ARG.LT.ARGMIN) GO TO 110
      SHALF=AK/ALOG(ARG)
      IF(ABS(SHALF-SHALF1).LT..0001) GO TO 102
      NIT=NIT+1
      IF(NIT.LT.11) GO TO 101
  102 S=SHALF**2
      SAV=.5*(S+SLOC)
      BP=1./(1.+BEE)
      GO TO 103
C     ------------------------------ LAMINAR FLOW
  110 AMRE=AM*RE
      FRE=EF*RE
      IF(ABS(AMRE).LT..01) GO TO 111
      AMRE=AMAX1(-60.,AMIN1(60.,AMRE))
      EXPMRE=EXP(AMRE)
      STORE=EXPMRE-1.-AMRE
      AMRESQ=AMRE*AMRE
      SRE=AMRE*(1.-STORE*FRE/AMRESQ)/(EXPMRE-1.)
      BP=SRE*STORE/AMRESQ+FRE*(STORE-.5*AMRESQ)/(AMRESQ*AMRE)
      GO TO 112
  111 SRE=(2.-FRE*(1.+AMRE*.33333))/(2.+AMRE)
      BP=SRE*(.5+AMRE*.16667)+FRE*(.16667+AMRE*.041667)
  112 IF(SRE.GT.TINY) GO TO 113
      SRE=TINY
      BP=.33333
  113 S=SRE/RE
      SAV=S
  103 T=S*RRUREF
C     -------------------- UNDER-RELAX BP
      BP=BPLAST+.5*(BP-BPLAST)
      BPLAST=BP
      OUT1=BP
      OUT2=T
      S1(KWALL)=SAV
      S2(KWALL)=RRUREF
      S3(KWALL)=UREF
      S4(KWALL)=AM
      S5(KWALL)=AMRE
      S6(KWALL)=RE
      GO TO 900
C-----------------------------------------------------------------------
CHAPTER C ------------------------ OTHER DEPENDENT VARIABLES
  200 SAV=S1(KWALL)
      RRUREF=S2(KWALL)
      UREF=S3(KWALL)
      AM=S4(KWALL)
      AMRE=S5(KWALL)
      RE=S6(KWALL)
      IF(MODEL.EQ.1) GO TO 210
C     -------------------------------- TURBULENT FLOW
      PRRAT=PRL(J)*RECPRT(J)
C     ----- THE FOLLOWING P-FUNCTION IS NOT TO BE USED FOR PRRAT.LT.0.5
      PJAY=9.*(PRRAT-1.)/PRRAT**.25
      S=SAV*RECPRT(J)/(1.+AMAX1(-.99999,PJAY*SQRT(ABS(SAV))))
      OUT2=S*RRUREF
      IF(J.NE.JH) GO TO 221
      OUT1=(H-1.)*.5*UREF**2
      GO TO 900
  221 OUT1=0.
      GO TO 900
C     ------------------------------ LAMINAR FLOW
  210 IF(ABS(AMRE).LT..01) GO TO 211
      S=AM/(EXP(PRL(J)*AMRE)-1.)
      GO TO 212
  211 S=RECPRL(J)/(RE+.5*RE*PRL(J)*AMRE)
  212 OUT2=S*RRUREF
      IF(J.NE.JH) GO TO 214
      OUT1=(PRL(JH)-1.)*.5*UREF**2
      GO TO 900
  214 OUT1=0.
C     ------------------------------------------------------------------
  900 IF(ITEST.EQ.1) GO TO 901
      WRITE(6,9000) J,I1,OUT1,OUT2
 9000 FORMAT(12H WALL TESTS,,3H J=,I3,4H I1=,I3,6H OUT1=,1PE10.3,
     1 6H OUT2=,E10.3)
  901 RETURN
      END

      SUBROUTINE PLOTS(X,IDIME,IMAX,XAXIS,Y,JDIME,JMAX,YAXIS)
      common/dbg/debug
      logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
C                                                                      *
c  subroutine for plotting j curves of y(i,j) against x(i).           *
c                                                                      *
c  x and y are scaled to the range 0. to 1., for plotting as           *
c    (y-ymin)/(ymax-ymin), the maximun and minimum values are printed  *
c    n.b. the x and y arrays must be redefined before each call plots. *
c  idime is the variable dimension for x.                              *
c  imax is the number of x values.                                     *
c  xaxis stores the name of the x-axis.                                *
c  jdime is the variable dimension for y.                              *
c  jmax is the number of curves to be plotted, (up to 30).             *
c  the array yaxis(j) stores the names of the curves,                  *
c    the fist character of each curve-name is used for plotting.       *
c  xsize alters the x-plot size by a factor of .2 to 1., in steps of .1*
c  ysize is the y-plot size factor fo .2 upwards in steps of .2        *
c    xsize=1., ysize=1. gives normal size plot.                        *
c                                                                      *
c***********************************************************************
      DIMENSION X(IDIME),Y(IDIME,JDIME),YAXIS(JDIME),
     1 A(101),YMAX(30),YMIN(30),DIGIT(11)
c!!!! equivalence deactivated. It causes YMIN to be over-written
cc      EQUIVALENCE (YMAX(1),A(1)),(YMIN(1),A(31))
      CHARACTER*1         DOT,CROSS,BLANK,DIGIT,A
      character*8 YAXIS
      save
      DATA DOT,CROSS,BLANK/1H.,1H+,1H /
     1,DIGIT/1h ,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H1/
      if(debug) write(6,*)'>>> plots'
C***** SET PLOT SIZE FACTOR
      XSIZE=0.6
      YSIZE=1.4
      ysize=0.2
C***** SCALING X-ARRAY TO RANGE 0 TO 100*XSIZE
      XR=100.*XSIZE
      XMAX=-1.E30
      XMIN=+1.E30
      IM=IMAX
      DO 1 I=1,IM
      XMAX=AMAX1(XMAX,X(I))
    1	XMIN=AMIN1(XMIN,X(I))
      S=XR/(XMAX-XMIN+1.E-30)
      DO 2 I=1,IM
    2 X(I)=(X(I)-XMIN)*S
c***** scaling y-array to range 0 to 50*ysize
      YR=50.*YSIZE
      JM=JMAX
      DO 4 J=1,JM
        YMAX(J)=-1.E30
        YMIN(J)=+1.E30
        DO 4 I=1,IM
cccc        write(6,*)'i,j,ymax(j),ymin(j)',i,j,ymax(j),ymin(j)
        YMIN(J)=AMIN1(YMIN(J),Y(I,J))
    4 YMAX(J)=AMAX1(YMAX(J),Y(I,J))
      DO 3 J=1,JM
cccc        write(6,*)'j,ymin(j)',j,ymin(j)
        S=YR/(YMAX(J)-YMIN(J)+1.E-30)
        DO 3 I=1,IM
    3 Y(I,J)=(Y(I,J)-YMIN(J))*S
c***** write curve names, with actual min and max values
      J=1
      L=IFIX(XR/10.)
      K=L
      GO TO 6
    5 J=J+L
      K=K+L
    6 K=MIN0(JM,K)
      WRITE(6,101) (YAXIS(I),I=J,K)
      WRITE(6,103) (YMIN(I),I=J,K)
      WRITE(6,104) (YMAX(I),I=J,K)
      IF(K-JM) 5,7,7
    7 WRITE(6,106)
c***** main loop - each pass produces a y-constant line
      IX=IFIX(XSIZE*10.)
      KX=IFIX(XR)+1
      IY=IFIX(YR*.1)
      KY=IFIX(YR)+1
      N=KY+1
      DO 40 M=1,KY
      L=N-M
      IF(L.EQ.1.OR.L.EQ.KY) GO TO 32
      GO TO 33
c***** put. or + along the x-axis
   32 DO 30 K=1,KX
   30 A(K)=DOT
      DO 31 K=1,KX,IX
   31 A(K)=CROSS
      GO TO 45
C***** PUT. OR + AND Y-VALUES ALONG Y-AXIS
   33 A(1)=DOT
      A(KX)=DOT
      K=L-1
   46 K=K-IY
      IF(K) 48,47,46
   47 A(1)=CROSS
      A(KX)=CROSS
   45 YL=FLOAT(L-1)/YR
      GO TO 35
   48 YL=-1.
C***** SEARCH FOR POINTS ALONG Y-CONSTANT LINE
   35 DO 42 J=1,JM
      DO 42 I=1,IM
      IF(IFIX(Y(I,J)+1.5).NE.L) GO TO 42
C***** POINT FOUND- ASSIGN SYMBOL- SEARCH FOR MORE
      NX=X(I)+1.5
      A(NX)=YAXIS(J)
   42 CONTINUE
C***** POINT Y-CONSTANT LINE
      IF(nint(yl).NE.-1) GO TO 37
      WRITE(6,106) (A(K),K=1,KX)
      GO TO 38
   37 WRITE(6,107) YL,(A(K),K=1,KX)
C***** FILL ARRAY A WITH BLANKS
   38 DO 49 K=1,KX
   49 A(K)=BLANK
   40 CONTINUE
C***** PRINT BLANK OR X-VALUE FOR X-AIXS
      L=1
      K=KX-1
      A(1)=DIGIT(1)
      DO 51 I=IX,K,IX
      L=L+1
      A(I)=DOT
   51 A(I+1)=DIGIT(L)
      A(K)=BLANK
      WRITE(6,106) (A(K),K=I,KX)
      WRITE(6,100) XAXIS,XMIN,XMAX
      RETURN
  100 FORMAT(1h ,12HABSCISSA IS ,A8,5H MIN=,1PE9.2,5H MAX=,E9.3)
  101 FORMAT(1h ,9HORDINATE ,12(A8,2X))
  103 FORMAT(1H ,7H   MIN ,1P11E10.2)
  104 FORMAT(1H ,7H   MAX ,1P11E10.2)
  106 FORMAT(6X,101A1)
  107 FORMAT(2H  ,F3.1,1X,101A1)
      END

      SUBROUTINE COMP
      common/dbg/debug
      logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
      DIMENSION A(20),B(20),C(20),CON(20),D(20),HCON(20),OMDIF(20)
      COMMON/COMA/
     1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
     2 EMU(20),F( 120),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
     3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
     4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
     5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
     6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
C
      EQUIVALENCE (A(1),DIF(1)),(C(1),SI(1)),(D(1),SIP(1))
      save
C
CHAPTER A --------------------------------------------------------------
C     INIT ------------- INIT ------------- INIT ------------- INIT
      if(debug) write(6,*)'>>> main'
      ENTRY INIT
      if(debug) write(6,*)'>>> init'
C     ------------------------- INITIAL VALUES AND DEFAULT VALUES
      NM1=N-1
      NM2=N-2
      NM3=N-3
      ISTEP=0
      IF(KRAD.EQ.3) NOVEL=1
      JUSTIN=ISTEP
      JUSTEX=ISTEP
      IFIN=1
      DXLAST=BIG
      DX=BIG
      PSII=0.
      BPE=1.
      BPI=1.
      Y(1)=0.
      DP=0.
      DO 13 I=1,N
      EMU(I)=0.
      CON(I)=0.
   13	R(I)=1.
      IF(NOVEL.NE.1) RETURN
      DO 16 I=1,N
   16	U(I)=1.
      RETURN
CCHPTER B --------------------------------------------------------------
C     GRID ------------- GRID ------------- GRID ------------- GRID 
      ENTRY GRID
      if(debug) write(6,*)'>>> grid'
      OMI=OM(2)
      OME=1.-OM(NM1)
      BOM(2)=.5*(OM(2)+OM(3))
      OMINT(1)=0.
      OMINT(2)=BOM(2)
      DO 202 I=3,NM2
      OMINT(I)=.5*(OM(I)+OM(I+1))
      BOM(I)=OMINT(I)-OMINT(I-1)
      OMDIF(I)=OM(I)-OM(I-1)
  202 HOMDFI=.5*OMDIF(3)
      BOM(NM1)=1.-OMINT(NM2)
      OMINT(NM1)=1.
      OMDIF(NM1)=OM(NM1)-OM(NM2)
      HOMDFE=.5*OMDIF(NM1)
      RETURN
CHAPTER C --------------------------------------------------------------
C     DISTAN ---------- DISTAN ---------- DISTAN ---------- DISTAN
      ENTRY DISTAN
      if(debug) write(6,*)'>>> distan'
      IF(NOVEL.NE.1) GO TO 220
      DO 224 I=1,N
  224	RECRU(I)=1./(RHO(I)+TINY)
      GO TO 222
  220	DO 221 I=1,N
      RECRU(I)=1./(RHO(I)*U(I)+TINY)
      IF(RECRU(I).GT.0.) GO TO 221
      IFIN=2
      WRITE(6,223) RECRU(I),I,ISTEP
  223	FORMAT(14H *** RECRU(I)=,1PE10.3,6H AT I=,I4,11H AND ISTEP=,I5,
     1 16H *** COMP DISTAN)
      stop
  221	CONTINUE
c     -------------------- calculation of y's and r's
  222	IF(KIN.EQ.1) GO TO 308
      RAT=RECRU(2)*RHO(1)*U(1)
      IF(KRAD.EQ.2) GO TO 307
      BPI=.33333+.66667*RAT
      GO TO 308
  307	BPI=(R(1)*(.83333*RAT+.16667)+.5*R(2)*(RAT+1.))/(R(1)+R(2))
  308	IF(KEX.EQ.1) GO TO 230
      RAT=RECRU(NM1)*RHO(N)*U(N)
      IF(KRAD.EQ.2) GO TO 327
      BPE=.33333+.66667*RAT
      GO TO 230
  327	BPE=(R(N)*(.83333*RAT+.16667)+.5*R(NM1)*(RAT+1.))/(R(N)+R(NM1))
c     ------------------- y's for plane flow
  230	STORE=OM(2)/BPI
      ADPEI(2)=(HOMDFI+STORE)*RECRU(2)
      Y(2)=PEI*RECRU(2)*STORE
      HPEI=.5*PEI
      DO 231 I=3,NM1
      ADPEI(I)=BOM(I)*RECRU(I)
  231	Y(I)=Y(I-1)+HPEI*OMDIF(I)*(RECRU(I)+RECRU(I-1))
      STORE=OME/BPE
      ADPEI(NM1)=(HOMDFE+STORE)*RECRU(NM1)
      Y(N)=Y(NM1)+PEI*RECRU(NM1)*STORE
C     ------------------------------------------------
      IF(KRAD-2) 270,240,280
C     ---------------- Y'S AND R'S FOR AXIAL SYMMETRY
  240 IF(CSALFA.LT.TINY) GO TO 260
C     ------------------------------ CSALFA NE 0.
      COSD2=.5*CSALFA
      TWDCOS=2./CSALFA
      IF(R(1).GT.TINY) GO TO 250
C     ------------------------------ R(1)=0.
      DO 242 I=2,N
      Y(I)=SQRT(ABS(Y(I)*TWDCOS))
  242 R(I)=Y(I)*CSALFA
      GO TO 270
C     ------------------------------ R(1) NE 0.
  250 R1D2=.5*R(1)
      R1D2SQ=R1D2*R1D2
      DO 251 I=2,N
      Y(I)=Y(I)/(R1D2+SQRT(ABS(R1D2SQ+COSD2*Y(I))))
  251 R(I)=R(1)+Y(I)*CSALFA
      GO TO 270
C     ------------------------------ CSALFA=0.
  260 RECR1=1./R(1)
      DO 261 I=2,N
      R(I)=R(1)
  261 Y(I)=Y(I)*RECR1
      GO TO 270
C     -------------------------------- POINT SYMMETRY, KRAD=3
  280 R1CUB=R(1)**3
      DO 281 I=2,N
      R(I)=(R1CUB+Y(I))**.3333333
  281 Y(I)=R(I)-R(1)
C     ------------------------------ GENERAL
  270 YI=Y(2)
      YE=Y(N)-Y(NM1)
      DO 273 I=1,NM1
  273 RECYDF(I)=1./(Y(I+1)-Y(I))
      IF(ITEST.EQ.1) RETURN
      WRITE(6,274) (RHO(I),I=1,N)
      WRITE(6,275) (RECRU(I),I=1,N)
      WRITE(6,276) (ADPEI(I),I=1,N)
      WRITE(6,277) (Y(I),I=1,N)
      WRITE(6,278) (R(I),I=1,N)
      WRITE(6,279) (RECYDF(I),I=1,N)
      RETURN
  274 FORMAT(18H0COMP DISTAN TESTS/8H RHO(I)=/(3X,1P6E11.3))
  275 FORMAT(10H RECRU(I)=/(3X,1P6E11.3))
  276 FORMAT(10H ADPEI(I)=/(3X,1P6E11.3))
  277 FORMAT(6H Y(I)=/(3X,1P6E11.3))
  278 FORMAT(6H R(I)=/(3X,1P6E11.3))
  279 FORMAT(11H RECYDF(I)=/(3X,1P6E11.3))
CHAPTER D --------------------------------------------------------------
C     SOLVE ---------- SOLVE ----------- SOLVE ---------- SOLVE
      ENTRY SOLVE
      if(debug) write(6,*)'>>> solve'
C     ---------------------------------------- PRELIMINARIES
      DXDPEI=DX/PEI
      CONST1=.5*DXDPEI
      ENT=ABS(RMI)+ABS(RME)
      IF(ENT.LE.TINY) GO TO 310
      HCONI=RMI*CONST1
      HCONDF=(RME-RMI)*CONST1
      DO 412 I=2,NM1
      HCON(I)=HCONI+HCONDF*OMINT(I)
  412 CON(I)=HCON(I)+HCON(I)
C     ---------------------------------------- COEFFICIENTS FOR U
  310 IF(NOVEL.EQ.1) GO TO 442
      J=0
C     ----- CALL SUBROUTINE PHYS AT ENTRY PHYSU
      CALL PHYSU
      IF(KRAD-2) 410,415,411
  410 DO 413 I=2,NM2
  413 DIFU(I)=CONST1*(EMU(I)+EMU(I+1))*RECYDF(I)
      GO TO 414
  415 CONST2=.5*CONST1
      DO 416 I=2,NM2
  416 DIFU(I)=CONST2*(R(I+1)+R(I))*(EMU(I)+EMU(I+1))*RECYDF(I)
      GO TO 414
  411 CONST3=.25*CONST1
      DO 419 I=2,NM2
  419 DIFU(I)=CONST3*(R(I+1)+R(I))**2*(EMU(I)+EMU(I+1))*RECYDF(I)
C     ----------------------------------------- A'S AND B'S
  414 IF(ENT.LE.TINY) GO TO 312
      DO 417 I=2,NM2
      A(I)=AMAX1(0.,DIFU(I)-HCON(I),-CON(I))
  417 B(I+1)=A(I)+CON(I)
      GO TO 314
  312 DO 315 I=2,NM2
      A(I)=DIFU(I)
  315 B(I+1)=A(I)
  314 TI=0.
      TE=0.
      IF(KIN.EQ.1) CALL WALL(1,BPI,TI)
      IF(KEX.EQ.1) CALL WALL(N,BPE,TE)
      B(2)=AMAX1((TI+RMI)*DXDPEI,0.)
      A(NM1)=AMAX1((TE-RME)*DXDPEI,0.)
cc      write(6,*)'TE,RME,DXDPEI',TE,RME,DXDPEI
C     ------------------------------------------ C'S AND D'S
      IF(MOMSOU.EQ.0) GO TO 431
      DO 418 I=2,NM1
      C(I)=U(I)*BOM(I)+SI(I)
  418 D(I)=A(I)+B(I)+BOM(I)  
      GO TO 432
  431 DO 433 I=2,NM1
      C(I)=U(I)*BOM(I)
  433 D(I)=A(I)+B(I)+BOM(I)
  432 CONTINUE
      IF(ITEST.EQ.1) GO TO 404
      WRITE(6,341) (DIFU(I),I=2,NM1)
      WRITE(6,342) (CON(I),I=2,NM1)
      WRITE(6,405) (A(I),I=2,NM1)
      WRITE(6,406) (B(I),I=2,NM1)
      WRITE(6,407) (C(I),I=2,NM1)
      WRITE(6,408) (D(I),I=2,NM1)
  341 FORMAT(23HOCOMP SOLVE TESTS FOR U/9H DIFU(I)=/(3X,1P6E11.3))
  342 FORMAT(8H CON(I)=/(3X,1P6E11.3))
  405 FORMAT(6H A(I)=/(3X,1P6E11.3))
  406 FORMAT(6H B(I)=/(3X,1P6E11.3))
  407 FORMAT(6H C(I)=/(3X,1P6E11.3))
  408 FORMAT(6H D(I)=/(3X,1P6E11.3))
  404 CONTINUE
C  ------------------------------- ADJUST FREE-BOUNDARY VALUES
      IF(KIN.EQ.2) U(1)=U(1)-DP*RECRU(1)
      IF(KEX.EQ.2) U(N)=U(N)-DP*RECRU(N)
C----------------------------   SOLVE FOR DOWNSREAM U'S
      C(2)=(B(2)*U(1)+C(2))/D(2)
      D(2)=A(2)/D(2)
      DO 421 I=3,NM1
      T=1./(D(I)-B(I)*D(I-1))
      D(I)=A(I)*T
  421 C(I)=(B(I)*C(I-1)+C(I))*T
      DO 422 IDASH=1,NM2
      I=N-IDASH
  422 U(I)=D(I)*U(I+1)+C(I)
      IF(KIN-2) 444,445,446
  444 TAUI=TI*U(2)/R(1)
      GO TO 445
  446 U(1)=U(2)
  445 IF(KEX-2) 447,440,448
  447 TAUE=TE*U(NM1)/R(N)
      GO TO 440
  448 U(N)=U(NM1)
  440 IF(ITEST.NE.1) WRITE(6,443) (U(I),I=1,N)
  443 FORMAT(6H U(I)=/(3X,1P6E11.3))
C     -------------------------------------------
C     ---------------------------------------------- F-SECTION
  442 IF(NF.LT.1) GO TO 481
      DO 480 J=1,NF
      IDJ=IDIMF*(J-1)
      I1J=1+IDJ
      I2J=2+IDJ
      INM1J=NM1+IDJ
      INJ=N+IDJ
C     ----- CALL SUBROUTINE PHYS AT ENTRY PHYSF
      CALL PHYSF
      TIF=0.
      FDIFI=0.
      TEF=0.
      FDIFE=0.
      IF(KIN.EQ.1) CALL WALL(1,FDIFI,TIF)
      IF(KEX.EQ.1) CALL WALL(N,FDIFE,TEF)
      IF(ITEST.EQ.1) GO TO 450
      WRITE(6,451) J,(DIF(I),I=2,NM1)
  451 FORMAT(24H COMP SOLVE TESTS FOR J=,I3/8H DIF(I)=/(3X,1P6E11.3))
C     ---------------------------------- COEFFICIENTS FOR F S
C     -------------------------------------------- A'S AND B'S
  450 IF(NEWPR.EQ.1) GO TO 337
      IF(ENT.LE.TINY) GO TO 335
      DO 484 I=2,NM2
      A(I)=AMAX1(0.,DIF(I)-HCON(I),-CON(I))
  484 B(I+1)=A(I)+CON(I)
      GO TO 337
  335 DO 338 I=2,NM2
      A(I)=DIF(I)
  338 B(I+1)=A(I)
  337 CONTINUE
      B(2)=AMAX1((TIF+RMI)*DXDPEI,0.)
      A(NM1)=AMAX1((TEF-RME)*DXDPEI,0.)
C     -------------------------------------------- C'S AND D'S
      GO TO (501,502,503), KSOURC
C     ------------------------------ KSOURC=1, GENERAL
  501 SI2=SI(2)
      SINM1=SI(NM1)
      DO 485 I=2,NM1
      IJ=I+IDJ
      D(I)=A(I)+B(I)+BOM(I)-SIP(I)
cc      write(6,*)'i,ij,f(ij)',i,ij,f(ij)
  485 C(I)=F(IJ)*BOM(I)+SI(I)
      GO TO 504
C     ------------------------------ KSOURC=2, NO SIP
  502 SI2=SI(2)
      SINM1=SI(NM1)
      DO 505 I=2,NM1
      IJ=I+IDJ
      D(I)=A(I)+B(I)+BOM(I)
cc      write(6,*)'i,ij,f(ij)',i,ij,f(ij)
  505 C(I)=F(IJ)*BOM(I)+SI(I)
      GO TO 504
C     ------------------------------ KSOURC=3, NO SIP OR SI
  503 SI2=0.
      SINM1=0.
      DO 506 I=2,NM1
      IJ=I+IDJ
      D(I)=A(I)+B(I)+BOM(I)
cc      write(6,*)'i,ij,f(ij)',i,ij,f(ij)
  506 C(I)=F(IJ)*BOM(I)
C     ------------------------------
  504 C(2)=C(2)-TIF*FDIFI*DXDPEI
      C(NM1)=C(NM1)-TEF*FDIFE*DXDPEI
      IF(KIN.GT.1) GO TO 486
      IF(IBIN(J).EQ.1) GO TO 486
      B(2)=0.
ccc      C(2)=F(120)*BOM(2)+SI2+RJTOTI(J)*DXDPEI
      C(2)=F(I2J)*BOM(2)+SI2+RJTOTI(J)*DXDPEI
      D(2)=D(2)-TIF*DXDPEI
  486 IF(KEX.GT.1) GO TO 491
      IF(IBEX(J).EQ.1) GO TO 491
      A(NM1)=0.
      C(NM1)=F(INM1J)*BOM(NM1)+SINM1-RJTOTE(J)*DXDPEI
      D(NM1)=D(NM1)-TEF*DXDPEI
  491 CONTINUE
      IF(ITEST.EQ.1) GO TO 464
      WRITE(6,405) (A(I),I=2,NM1)
      WRITE(6,406) (B(I),I=2,NM1)
      WRITE(6,407) (C(I),I=2,NM1)
      WRITE(6,408) (C(I),I=2,NM1)
C     ------------------------------ SOLVE FOR DOWNSTREAM F'S
  464 C(2)=(B(2)*F(I1J)+C(2))/D(2)
      D(2)=A(2)/D(2)
      DO 465 I=3,NM1
      T=1./(D(I)-B(I)*D(I-1))
      D(I)=A(I)*T
  465 C(I)=(B(I)*C(I-1)+C(I))*T
      DO 466 IDASH=1,NM2
      I=N-IDASH
      IJ=I+IDJ
  466 F(IJ)=D(I)*F(IJ+1)+C(I)
C     --------------------------- ADJUST F(1,J) AND F(N,J)
      IF(KIN-2) 467,460,469
  467 IF(IBIN(J).EQ.1) GO TO 468
      F(I1J)=FDIFI+F(I2J)+(RJTOTI(J)-F(I1J)*RMI)/TIF
      GO TO 460
  469 F(I1J)=F(I2J)
      GO TO 460
  468 RJTOTI(J)=TIF*(F(I1J)-F(I2J)-FDIFI)+RMI*F(I1J)
  460 IF(KEX-2) 471,470,473
  471 IF(IBEX(J).EQ.1) GO TO 472
      F(INJ)=FDIFE+F(INM1J)-(RJTOTE(J)-F(INJ)*RME)/TEF
      GO TO 470
  473 F(INJ)=F(INM1J)
      GO TO 470 
  472 RJTOTE(J)=TEF*(F(INM1J)+FDIFE-F(INJ))+RME*F(INJ)
  470 IF(ITEST.EQ.1) GO TO 480
      WRITE(6,476) J,(F(I+IDJ),I=1,N)
  476 FORMAT(6H  F(I,,I2,1H)/(3X,1P6E11.3))
  480 CONTINUE
C     ------------------------------------------
  481 XU=XD
      PSII=PSII-RMI*DX
      PSIE=PSIE-RME*DX
      PEI=PSIE-PSII
      ISTEP=ISTEP+1
      END
      SUBROUTINE PHYS
      common/dbg/debug
      logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
      COMMON/COMA/
     1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
     2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
     3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
     4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
     5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
     6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
       COMMON/COMB/
     1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
     2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
     3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
     4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
     5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
     6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
     7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
     8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
     9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
     1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
     2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
C
      DIMENSION DUDY(20),EL(20),YEDGE(6)
      save
C
c#### should be in blockdata
ccccc      DATA KUDIF/-1/
C     ------------------------------------------------------------------
CHAPTER A ---------- PHYSU ---------- PHYSU ---------- PHYSU -----------
      if(debug) write(6,*)'>>> phys'
      ENTRY PHYSU
      if(debug) write(6,*)'>>> physu'
C     ---------------------------------------- LAMINAR VISCOSITY
C     SQUARE-ROOT FORMULA, WITH WEIGHTING ACCORDING TO MASS FRACTION
      DO 110 I=1,N
  110	EMU(I)=(VISFU*F(I,JF)+VISOX*F(I,JOX)+VISPR*F(I,JPR))*
     1                                               SQRT(F(I,JTE))
      IF(MODEL.EQ.1) GO TO 209
C     ------------------------------------------------------------------
C     ----------------------------- MIXING LENGTH MODEL OF TURBULENCE
      IF(KUDIF.EQ.ISTEP) GO TO 102
      UMAX=U(1)
      UMIN=U(1)
      DO 101 I=2,N
      UMAX=AMAX1(UMAX,U(I))
  101	UMIN=AMIN1(UMIN,U(I))
      UDIF=UMAX-UMIN
  102	HUDIF=.5*UDIF
      DUDYMN=FR*UDIF/Y(N)
C
      DO 105 I=2,NM1
  105	DUDY(I)=ABS(U(I+1)-U(I-1))/(Y(I+1)-Y(I-1))
      K=1
      EX=DUDY(2)-DUDYMN
      IF(EX.LT.0.) GO TO 103
      YEDGE(K)=0.
      K=2
  103	DO 104 I=3,NM1
      EXL=EX
      EX=DUDY(I)-DUDYMN
      IF(EX*EXL.GE.0.) GO TO 104
      YEDGE(K)=.5*(Y(I)+Y(I-1))
      IF(K.EQ.6) GO TO 107
      K=K+1
  104	CONTINUE
      IF(EX.LT.0.) GO TO 108
      YEDGE(K)=Y(N)
      IF(K.EQ.6) GO TO 107
      K=K+1
  108	CONTINUE
      DO 106 KAY=K,6
  106	YEDGE(KAY)=Y(N)
  107	CONTINUE
      EL12=(YEDGE(2)-YEDGE(1))*ALMG
      EL34=(YEDGE(4)-YEDGE(3))*ALMG
      EL56=(YEDGE(6)-YEDGE(5))*ALMG
      EL23=.5*(EL12+EL34)
      EL45=.5*(EL34+EL56)
      ASSIGN 119 TO K
      DO 130 I=2,NM1
      YVALUE=Y(I)
      GO TO K, (119,121,123,125,127,129)
  119	IF(YVALUE.LT.YEDGE(1)) GO TO 120
      ASSIGN 121 TO K
  121	IF(YVALUE.LT.YEDGE(2)) GO TO 122
      ASSIGN 123 TO K
  123	IF(YVALUE.LT.YEDGE(3)) GO TO 124
      ASSIGN 125 TO K
  125	IF(YVALUE.LT.YEDGE(4)) GO TO 126
      ASSIGN 127 TO K
  127	IF(YVALUE.LT.YEDGE(5)) GO TO 128
      ASSIGN 129 TO K
      GO TO 129
  120	EL(I)=0.
      GO TO 130
  122	EL(I)=EL12
      GO TO 130
  124	EL(I)=EL23
      GO TO 130
  126	EL(I)=EL34
      GO TO 130
  128	EL(I)=EL45
      GO TO 130
  129	EL(I)=EL56
C     ------------------------------ UPPER LIMITS TO MIXING LENGTH
  130	EL(I)=AMIN1(EL(I),HUDIF/(DUDY(I)+TINY))
      IF(KIN.NE.1) GO TO 141
      DO 142 I=2,NM1
  142	EL(I)=AMIN1(EL(I),AK*Y(I))
  141	IF(KEX.NE.1) GO TO 143
      DO 144 I=2,NM1
  144	EL(I)=AMIN1(EL(I),AK*(Y(N)-Y(I)))
  143	CONTINUE
C
C     ------------------------------------------- VISCOSTITES
C  -------------------------------------------TURBULENT CONTRIBUTION
cccc  200   DO 201 I=2,NM1
      DO 201 I=2,NM1
      DUDYL=DUDY(I)*EL(I)
      UDMIN=UFAC*U(I)
      DUDYL=AMAX1(DUDYL,UDMIN)
      EMUT=RHO(I)*EL(I)*DUDYL
C     ----- SIMPLE ADDITION OF THE TURBULENT AND LAMINAR CONTRIBUTIONS
      EMU(I)=EMU(I)+EMUT
  201	CONTINUE
C
C     ------------------------------------------------MOMENTUM SOURCE
  209	AGRVDX=AGRAV*DX
      RPRLST=1.
      MOMSOU=1
      IF(ABS(DP).GT.TINY) GO TO 204
      IF(ABS(AGRAV).GT.TINY) GO TO 204
      MOMSOU=0
      RETURN
  204 continue
      DO 210 I=2,NM1
  210	SI(I)=ADPEI(I)*(AGRVDX*(RHO(N)-RHO(I))-DP)
      GO TO 9000
C
C     ------------------------------------------------------------------
chapter b ---------- physf ---------- physf ---------- physf ----------- 
      ENTRY PHYSF
      if(debug) write(6,*)'>>> physf'
      IF(MODEL.EQ.2) GO TO 312
      RECPR=RECPRL(J)
      GO TO 310
  312	RECPR=RECPRT(J)
  310	NEWPR=1
      IF(ABS(RECPR-RPRLST).LT.1.E-10) GO TO 314
      NEWPR=2
      DO 313 I=2,NM2
  313	DIF(I)=DIFU(I)*RECPR
      RPRLST=RECPR
c     ------------------------------------------ kinetic heating source
  314	IF(J.NE.JH) GO TO 3000
      IF(ABS(RECPR-1.).LT.1.E-10) GO TO 320
      IF(NOVEL.NE.1) GO TO 321
  320 KSOURC=3
      RETURN
  321 DUSQP=0.
      USQP=U(2)**2
      DO 322 I=2,NM2
        USQ=U(I+1)**2
        DUSQ=(DIFU(I)-DIF(I))*(USQ-USQP)
        SI(I)=.5*(DUSQ-DUSQP)
        si(i)=0.0
        DUSQP=DUSQ
  322 USQP=USQ
      KSOURC=2
      SI(NM1)=-.5*DUSQP
      si(nm1)=0.0
      GO TO 9000
C
C     ------------------------------------------------ FUEL SOURCE
 3000	IF(J.NE.JF) GO TO 4000
      IF(INERT.EQ.2) GO TO 342
      KSOURC=3
      RETURN
  342	KSOURC=1
      IF(MODEL.NE.1) GO TO 352
      T1=DX*PREEXP*PRESS**2
      T2=.5/STOICH
      DO 344 I=2,NM1
      FUBRNT=T2*(ABS(F(I,JP))-F(I,JP))
      IF(F(I,JF).GT.FUBRNT) GO TO 346
      SIP(I)=-BIG
      GO TO 344
  346	EXPO=EXP(-ARRCON/F(I,JTE))
      F(I,JOX)=AMAX1(0.,F(I,JP)+F(I,JF)*STOICH)
      TERM=-T1*EXPO*ADPEI(I)*F(I,JOX)
      SIP(I)=TERM/(1.-FUBRNT/F(I,JF))
  344	SI(I)=-SIP(I)*FUBRNT
      GO TO 9000
  352	T2=.5/STOICH
      T3=1./(PHIB-PHIC)
      T4=-PHIC*T3
C     ----- RATE CONTROLLED BY EDDY BREAK-UP
      CEBUDX=CEBU*DX
      DO 353 I=2,NM1
      FUBRNT=T2*(ABS(F(I,JP))-F(I,JP))
      FUUNBT=T3*F(I,JP)+T4
      IF(F(I,JF).LT.FUUNBT) then
        SIP(I)=-ADPEI(I)*CEBUDX*DUDY(I)*RHO(I)*(FUUNBT-F(I,JF))/
     1                                             (FUUNBT-FUBRNT+TINY)
      else
        SIP(I)=-BIG
      endif  
      SI(I)=-SIP(I)*FUBRNT
cc      write(6,*)'i,fubrnt,fuunbt',i,fubrnt,fuunbt
cc      write(6,*)'i,si(i),sip(i)',i,si(i),sip(i)
  353 continue
      GO TO 9000
C
C     ---------------------------------------------- JP (PHI)
 4000	IF(J.NE.JP) GO TO 9000
      KSOURC=3
      RETURN
C
 9000	IF(ITEST.EQ.1) RETURN
      WRITE(6,9001) J,(SI(I),I=2,NM1)
      WRITE(6,9002) (SIP(I),I=2,NM1)
 9001	FORMAT(18H PHYS TESTS FOR J=,I3/7H SI(I)=/(3X,1P6E11.3))
 9002	FORMAT(8H SIP(I)=/(3X,1P6E11.3))
      END
c-----------------------------------------------------------------------
      SUBROUTINE OUTPUT
      common/dbg/debug
      logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
      COMMON/COMA/
     1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
     2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
     3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
     4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
     5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
     6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
       COMMON/COMB/
     1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
     2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
     3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
     4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
     5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
     6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
     7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
     8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
     9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
     1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
     2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
C   
       DIMENSION LAB(6),OUT(6),TITLE(3,4),
     1 XLPLOT(100),YLAXIS(12),YLPLOT(100,12),
     2 XTPLOT(20),YTAXIS(4),YTPLOT(20,4),
     3 DFE(3),DFI(3),FLUX(3),STANE(3),STANI(3)
      character*10 lab
      character*8 ylaxis,ytaxis,title
      character*4 xtaxis,xlaxis
      save
C
chapter a ----------------------- initial data for printout ------------
c     -------------------- cross-stream output (profile) data ----------
c     ----- assign kout= no. of variables, and output labels lab(k)
      DATA KOUT/6/
      DATA (LAB(K),K=1,6)/'Y','U VEL','TEMP','FUEL','OXYG','OX-FU.STOI'/
C
c     -------------------- transverse (cross-stream) plot data
c     ----- assign nyt= no. of variables to be plotted
c     --- insert dimensions, ensure that itdim.ge.n.and.jtdim.ge.nyt
      DATA NYT/4/, ITDIM,JTDIM/20,4/
c     ----- assign labels for plot axes
      DATA XTAXIS/'Y(I)'/
      DATA (YTAXIS(K),K=1,4)/'U VEL   ','TEMP    ','FUEL    ',
     1                       'OXYG    '/
C
c     -------------------- longitudinal (down-stream) plot data
c     ----- assign nyl= no. of variables to be plotted
c     -- insert dimensions, ensure that ildim.ge.lastep.and.jldim.ge.nyl
      DATA NYL/12/ILDIM,JLDIM/100,12/
c!!!! two lines were missing
c     ----- assign labels for plot axes
      DATA XLAXIS/'XU'/
      DATA (YLAXIS(K),K=1,12)/'U(1)    ','T(1)    ','FU(1)   ',
     1                        'OX(1)   ','N,R OR Y','1,R(1)  ',
     1                        '2,PEI   ','3,RME   ','4,FLUXFU',
     1                        '5,DPDX  ','6,RATE  ','7,FACE  '/
C
C     -------------------- TITLE DATA
      DATA TITLE/'AXI-SYMM','ETRICAL ','FLOW     ',
     1 'PLANE FL','UW      ',' ',
     2 'RADIALLT','-OUTWARD',' FLOW   ',
     3 'VARIABLE',' CSALFA ',' '/
C
      if(debug) write(6,*)'>>> output'
CHAPTER B ----------------------------- HEADINGS -----------------------
      IF(ISTEP.GT.0) GO TO 1102
      WRITE(6,1103) (TITLE(I,KIND),I=1,3)
 1103   FORMAT(1h ,23HGENMIX, FEBRUARY 1977, ,3A8/1H ,
     1 'COMBUSTION OF METHANE AND AIR IN A DUCT.')
C
      PRESS1=PRESS
      TEM=.5*(R(1)+R(N))
      EMU1=(VISFU*F(1,JF)+VISOX*F(1,JOX)+VISPR*F(1,JPR))*
     1                                                 SQRT(F(1,JTE))
C
      REY=PEI/(EMU1*TEM)
      EQRAT=0.0
      IF(INERT.NE.1) EQRAT=FLOB*STOICH/(FLOC+TINY)/OXC
      AMACH=SQRT(PEI*UB/(GAMMA*PRESS*TEM))
C
      WRITE(6,1013) KASE,IRUN,KIND,KRAD,CSALFA,MODEL,INERT,NOVEL
 1013   FORMAT(1h ,5H KASE,5H IRUN,5H KIND,5H KRAD,7H CSALFA,6H MODEL,
     1 6H INERT,6H NOVEL/1X,I4,3I5,F7.3,3I6/)
C
      WRITE(6,1014) OMPOW, (OM(I),I=1,N)
 1014   FORMAT(1h ,18H OM(I), FOR OMPOW=,F6.3/(1X,1P6E11.3))
C
      WRITE(6,1010)
     1 HEX0,XHEX0,AHEX,BHEX,CHEX,
     2 HIN0,XHIN0,AHIN,BHIN,CHIN,
     3 UEX0,XUEX0,AUEX,BUEX,CUEX,
     4 XEND,XOUT,XULAST,HDIV,AGRAV
 1010   FORMAT(/1h ,
     1 4X,4HHEX0,6X,5HXHEX0,7X,4HAHEX,7X,4HBHEX,7X,4HCHEX/1X,1P5E11.3/
     2 5X,4HHIN0,6X,5HXHIN0,7X,4HAHIN,7X,4HBHIN,7X,4HCHIN/1X,1P5E11.3/
     3 5X,4HUEX0,6X,5HXUEX0,7X,4HAUEX,7X,4HBUEX,7X,4HCUEX/1X,1P5E11.3/
     4 5X,4HXEND,7X,4HXOUT,5X,6HXULAST,7X,4HHDIV,6X,5HAGRAV/1X,1P5E11.3)
C
      WRITE(6,1011) UA,UB,UC,UD,TA,TB,TC,TD,
     2 PRESS,PREEXP,REY,EQRAT,AMACH,ULIM,PEILIM
 1011   FORMAT(1h ,4X,2HUA,7X,2HUB,7X,2HUC,7X,2HUD,
     1  7X,2HTA,7X,2HTB,7X,2HTC,7X,2HTD/1X,8F9.3/
     2 4X,5HPRESS,3X,6HPREEXP,6X,3HREY,4X,5HEQRAT,4X,5HAMACH,5X,4HULIM,
     3 3X,6HPEILIM/1X,1P7E9.2)
C
chapter c ---------------- compute output required at each step --------
 1102	CONTINUE
      UBAR=0.
      DO 110 I=2,NM1
  110	UBAR=UBAR+BOM(I)*U(I)
      UFLUX=PEI*UBAR
      DO 115 J=1,NF
      FLUX(J)=0.
      DO 116 I=2,NM1
  116	FLUX(J)=FLUX(J)+BOM(I)*F(I,J)
  115	FLUX(J)=PEI*FLUX(J)
C
      DO 117 J=1,NF
      DFI(J)=FLUX(J)/PEI-F(I,J)
  117	DFE(J)=DFI(J)+F(1,J)-F(N,J)
      UFLUX=UFLUX-PSIE*U(N)+U(1)*PSII
      FLUX(JH)=FLUX(JH)-PSIE*ENTHD+PSII*ENTHA
      FLUX(JP)=FLUX(JP)-PSIE*PHID+PSII*PHIA
      FLUX(JF)=FLUX(JF)-PSIE*FUD+PSII*FUA
      PRESSD=PRESS/PRESS1-1.
C
      IF(ISTEP.EQ.0.OR.ILPLOT.EQ.1) GO TO 1105
c     ----- assign values for downstream plot
      XLPLOT(ISTEP)=XU
      YLPLOT(ISTEP,1)=U(1)
      YLPLOT(ISTEP,2)=F(1,JTE)    
      YLPLOT(ISTEP,3)=F(1,JF)
      YLPLOT(ISTEP,4)=F(1,JOX)
      IF(KIND-1) 111,111,114
  111	YLPLOT(ISTEP,5)=R(N)
      GO TO 113
  114	YLPLOT(ISTEP,5)=Y(N)
  113	CONTINUE
      YLPLOT(ISTEP,6)=R(1)
      YLPLOT(ISTEP,7)=PEI
      YLPLOT(ISTEP,8)=RME
      YLPLOT(ISTEP,9)=FLUX(JF)
      YLPLOT(ISTEP,10)=DPDX
      YLPLOT(ISTEP,11)=RATE
      YLPLOT(ISTEP,12)=FACE
 1105	CONTINUE
C
c     ----- tests for printout
c     ----- iprint=1 gives single (station) variables,
c           iprint=2 adds the array (profile) variables,
c           iprint-3 adds the cross-stream plots.
      IPRINT=0
      IF(MOD(ISTEP,NSTAT).EQ.0) IPRINT=1
      IF(MOD(ISTEP,NPROF).EQ.0) IPRINT=2
      IF(ISTEP.EQ.0) GO TO 1020
      IF(MOD(ISTEP,NPLOT).EQ.0
     1 .OR.ISTEP.EQ.JUSTEX.OR.ISTEP.EQ.JUSTIN
     2 .OR.ITEST.NE.1.OR.IFIN.NE.1) IPRINT=3
 1020	IF(IPRINT.EQ.0) RETURN
C
chapter d ------------------------------ station variables -------------
      call writbl
      WRITE(6,1030) XU,ISTEP,
     1 JUSTIN,JUSTEX,DX,PRESSD,
     2 KIN,KEX,DXY,DPDX,
     3 PSII,PSIE,DXRE,PEI,
     4 RMI,RME,DXINC,
     5 R(1),R(N),DXPSI,
     5 (J,J=1,4),UFLUX,(FLUX(J),J=1,NF)
 1030   FORMAT(1h ,5H***  ,3HXU=,1PE10.3,2X,6HISTEP=,I5/
     1 2X,7HJUSTIN=,I10,1X,7HJUSTEX=,I10,5X,3HDX=,1PE10.3,
     1 8H PRESSD=,E10.3/
     2 5X,4HKIN=,I10,4X,4HKEX=,I10,4X,4HDXY=,E10.3,3X,5HDPDX=,E10.3/
     3 4X,5HPSII=,E10.3,3X,5HPSIE=,E10.3,3X,5HDXRE=,E10.3,
     3 4X,4HPEI=,E10.3/
     4 5X,4HRMI=,E10.3,4X,4HRME=,E10.3,2X,6HDXINC=,E10.3/
     5 4X,5HR(1)=,E10.3,3X,5HR(N)=,E10.3,2X,6HDXPSI=,E10.3/
     5 23X,3HJ =,4(I6,5X)/7H UFLUX=,E10.3,9H FLUX(J)=,(4E11.3))
C
      IF(ISTEP.EQ.0) GO TO 1042
      UREF=UBAR
      RUREF=PEI/((R(I)+R(N))*.5*Y(N))
      URUREF=1./(UREF*RUREF)
C
      IF(KIN-2) 1061,1062,1063
 1061	TAUID=TAUI*URUREF
      DO 1025 J=1,NF
 1025	STANI(J)=(RJTOTI(J)-F(1,J)*RMI)/(R(1)*DFI(J)*RUREF+TINY)
      WRITE(6,1029) TAUID,(STANI(J),J=1,NF)
 1029	FORMAT(1H ,6HTAUID=,1PE10.3,10H STANI(J)=,(4E11.3))
      GOTO 1063
 1062	WRITE(6,1069) FACI,RATI
 1069	FORMAT(1H ,6H FACI=,1PE10.3,6H RATI=,E10.3)
 1063	IF(KEX-2) 1081,1082,1044
cccc 1081   TAUED=TAUED*URUREF
 1081   TAUED=TAUE*URUREF
      DO 1027 J=1,NF
 1027   STANE(J)=(RJTOTE(J)-F(N,J)*RME)/(R(N)*DFE(J)*RUREF+TINY)
      WRITE(6,1028) TAUED,(STANE(J),J=1,NF)
 1028	FORMAT(1H ,6HTAUED=,1PE10.3,10H STANE(J)=,(4E11.3))
      GO TO 1044
 1082	WRITE(6,1089) FACE,RATE
 1089	FORMAT(1H ,6H FACE=,1PE10.3,5H DA2=,E10.3)
      GO TO 1042
 1044	WRITE(6,1047) DA1,DA2
 1047	FORMAT(5H DA1=,1PE10.3,5H DA2=,E10.3)
C
CCHAPTER E ---------------------------- CROSS-STREAM PROFILES ----------
 1042 IF(IPRINT.EQ.1) GO TO 1050
C
      WRITE(6,1099) (LAB(K),K=1,KOUT)
      DO 1091 I=1,N
      OUT(1)=Y(I)
      OUT(2)=U(I)
      OUT(3)=F(I,JTE)
      OUT(4)=F(I,JF)
      OUT(5)=F(I,JOX)
      OUT(6)=F(I,JP)
c     ------------------------------ print profiles
 1091 WRITE(6,1098) I,(OUT(K),K=1,KOUT)
C
      IF(IPRINT.LT.3.OR.ITPLOT.EQ.1) GO TO 1050
c     ----- assign cross-stream plots
      DO 1073 I=1,N
      XTPLOT(I)=Y(I)
      YTPLOT(I,1)=U(I)
      YTPLOT(I,2)=F(I,JTE)
      YTPLOT(I,3)=F(I,JF)
      YTPLOT(I,4)=F(I,JOX)
 1073	CONTINUE
c     --------------------- cross-stream plot output
      call writbl
      WRITE(6,1096) XU,ISTEP
 1096   FORMAT(19H CROSS-STREAM PLOT,,4H XU=,1PE10.3,7H ISTEP=,I4)
      CALL PLOTS(XTPLOT,ITDIM,N,XTAXIS,YTPLOT,JTDIM,NYT,YTAXIS)
C
chapter f -------------------------- return or terminate ---------------
 1050	IF(IFIN.EQ.1) RETURN
C
      call writbl
      WRITE(6,112) ISTEP,LASTEP,XU,XULAST,IFIN
  112	FORMAT(14H0TERMINATED AT/7H ISTEP=,I5,8H LASTEP=,I5,
     1 4H XU=,1PE11.3,8H XULAST=,E11.3,6H IFIN=,I3)
      IF (ILPLOT.EQ.1) RETURN
      call writbl
c     --------------------- downstream plot output
      WRITE(6,1054) XU,ISTEP
 1054   FORMAT(18H DOWN-STREAM PLOT,,4H XU=,1PE10.3,7H ISTEP=,I4)
      CALL PLOTS(XLPLOT,ILDIM,ISTEP,XLAXIS,YLPLOT,JLDIM,NYL,YLAXIS)
 1098	FORMAT(1H ,I3,1P10E11.3)
 1099   FORMAT(1h ,2X,2HI ,10A11)
      END
      subroutine writbl
      write(6,*)' '
      end
      
c