c

C.... File name ................. GXDENS.HTM .................... 021110
C
      FUNCTION DENSITY(I)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/prpcmn'
      COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(39),
     1            ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
     1      /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
      COMMON/HBASE/IH01,IH02,KH01,KH01H,KH01L,KH02,KH02H,KH02L,L0H012
      common/dbe/dbear
      logical dbear
C
c
      IF(IGRND.EQ.-1) THEN ! default when RHO1 or 2 is greater than zero
        DENSITY=PRPRTY
      ELSEIF(IGRND.LE.2) THEN ! when RHO1 or 2=GRND1 or GRND2
c
C.... Density is a linear function or the reciprocal of a linear function
C     of either a scalar (enthalpy or other scalar); or two scalars; or
C     three:
c
        IF(L0SCAL.GT.0) DENSITY=RHOAG   + RHOBG*F(L0SCAL+I)
        IF(L0CA.GT.0)   DENSITY=DENSITY + RHOCG*F(L0CA+I)
        IF(L0CB.GT.0)   DENSITY=DENSITY + RHODG*F(L0CB+I)
c
        IF(IGRND.EQ.2) DENSITY=1./DENSITY
c
      ELSEIF(IGRND.EQ.3) THEN ! when RHO1 or 2=GRND1
c
C.... Density for the compressible flow and also shallow-water theory:
cccc        call writ3r('rhoag  ',rhoag,'rhobg   ',rhobg,'rhocg   ',rhocg)
        DENSITY=RHOAG * (AMAX1(0.,ABSPRES(I)))**RHOBG + RHOCG
c
C.... Multiplication by C3 is useful when two compressible gases are
C     present in the same domain. C3 has of course to be computed so as
C     to "travel" with the gas which it marks:
c
        IF(L0SCAL.NE.0) DENSITY=DENSITY*F(L0SCAL+I)
c
C.... For shallow-water theory with NZ > 1, and compressibility in the
C     upper (lower-IZ) cell only, the density is set to -RHO2 except in
C     topmost cells (IZ=1):
c
        IF(SHALWT) DENSITY=- RHO2
c
      ELSEIF(IGRND.EQ.4) THEN
c
C.... Density is a linear function of temperature:
c
        DENSITY=RHOAG + RHOBG*(F(L0TEM+I)+TEMP0)
c
      ELSEIF(IGRND.EQ.5) THEN
c
C.... Density according to the ideal-gas law:
        ABST=ABSTEMP(I)                         ! get absolute temperature
        IF(L0SCAL.GT.0.AND.RHOAG.GT.0.0) THEN   ! for a mixture of 2 gases
          SPEVOL=RHOAG + (RHOBG-RHOAG) * F(L0SCAL+I)  ! the scalar is C1 or
          DENSITY=ABSPRES(I) / (SPEVOL * ABST)         ! C2, according to phase
        ELSE                                            ! RHOAG and RHOBG are
          DENSITY=RHOBG * ABSPRES(I) / ABST            ! the ideal-gas constant
        ENDIF                                           ! of the gases
c
       ELSEIF(IGRND.EQ.6) THEN
c
C.... Ideal-gas law for mixture of the three chemical 'species' fuel,
C     oxidant and product involved in the simple chemical reaction
C     model (Parameters RHOA,RHOB and RHOC are the molecular weights
C     of fuel, oxidant and products respectively):
c
        ABST=ABSTEMP(I)
        IF(ABST.GT.0.0) THEN
          RMIX=(F(L0FUEL+I)/RHOAG + F(L0OXID+I)/RHOBG +
     1                                    F(L0PROD+I)/RHOCG) * GASCON
          DENSITY=ABSPRES(I)/(RMIX * ABST)
        ELSE
          call writ40('absolute temperature less than zero')
          call writ40('density set to 1.0                 ')
          density=1.0
        ENDIF
c
      ELSEIF(IGRND.EQ.7) THEN
c
C.... Ideal-gas law for a mixture of the 7 gases involved in a simple
C     model of the combustion of hydrocarbons in air.
Click here for explanation
C.... Mass fraction of anything but air is FNOTAIR.
C!!!! Note that prior to the generalization of March 2002, it has been
c     that any material injected which did not have a finite MIXF value
c     was air, with .232 O2 content and .768 N2 content.
c     This presumption is preserved here.
c
C.... FO, FC and FH are the mass fractions of the elements O, C and H,
C     regardless of which compound they are contained in.
        FC=0.0 ; FH=0.0 ; FO=0.0 ; FN=0.0  ! initialise to zero
        FNOTAIR=0.0
        IF(L0MIXF.NE.0) THEN
          FMXF0=F(L0MIXF+I)        ! mass fraction of material from inlet 0
          FH=FH + RHOBG * FMXF0       ! stream * hydrogen content (rhobg)
          FC=FC + RHOAG * FMXF0       !          carbon content   (rhoag)
          FN=FN + RHOCG * FMXF0       !          nitrogen content (rhocg)
                                       ! inlet 0 contains no oxygen
          FNOTAIR=FNOTAIR + FMXF0
        ENDIF
        IF(L0MXF1.NE.0) THEN
          FMXF1=F(L0MXF1+I)        ! mass fraction of material from inlet 1
          FH=FH + FHMXF1 * FMXF1      ! stream * hydrogen content (FH1)
          FC=FC + FCMXF1 * FMXF1      !          carbon content   (FC1)
          FO=FO + FOMXF1 * FMXF1      !          oxygen content   (FO1)
          FN=FN + FNMXF1 * FMXF1      !          nitrogen content (FN1)
          FNOTAIR=FNOTAIR + FMXF1
        ENDIF
        IF(L0MXF2.NE.0) THEN
          FMXF2=F(L0MXF2+I)        ! mass fraction of material from inlet 2
          FH=FH + FHMXF2 * FMXF2      ! stream * hydrogen content (FH2)
          FC=FC + FCMXF2 * FMXF2      !          carbon content   (FC2)
          FO=FO + FOMXF2 * FMXF2      !          oxygen content   (FO2)
          FN=FN + FNMXF2 * FMXF2      !          nitrogen content (FN2)
          FNOTAIR=FNOTAIR + FMXF2
        ENDIF
        IF(L0MXF3.NE.0) THEN
          FMXF3=F(L0MXF3+I)        ! mass fraction of material from inlet 3
          FH=FH + FHMXF3 * FMXF3      ! stream * hydrogen content (FH3)
          FC=FC + FCMXF3 * FMXF3      !          carbon content   (FC3)
          FO=FO + FOMXF3 * FMXF3      !          oxygen content   (FO3)
          FN=FN + FNMXF3 * FMXF3      !          nitrogen content (FN3)
          FNOTAIR=FNOTAIR + FMXF3
        ENDIF
        IF(L0MXF4.NE.0) THEN
          FMXF4=F(L0MXF4+I)        ! mass fraction of material from inlet 4
          FH=FH + FHMXF4 * FMXF4      ! stream * hydrogen content (FH4)
          FC=FC + FCMXF4 * FMXF4      !          carbon content   (FC4)
          FO=FO + FOMXF4 * FMXF4      !          oxygen content   (FO4)
          FN=FN + FNMXF4 * FMXF4      !          nitrogen content (FN4)
          FNOTAIR=FNOTAIR + FMXF4
        ENDIF
        IF(L0MXF5.NE.0) THEN
          FMXF5=F(L0MXF5+I)        ! mass fraction of material from inlet 5
          FH=FH + FHMXF5 * FMXF5      ! stream * hydrogen content (FH5)
          FC=FC + FCMXF5 * FMXF5      !          carbon content   (FC5)
          FO=FO + FOMXF5 * FMXF5      !          oxygen content   (FO5)
          FN=FN + FNMXF5 * FMXF5      !          nitrogen content (FN5)
          FNOTAIR=FNOTAIR + FMXF5
        ENDIF
        IF(L0MXF6.NE.0) THEN
          FMXF6=F(L0MXF6+I)        ! mass fraction of material from inlet 6
          FH=FH + FHMXF6 * FMXF6      ! stream * hydrogen content (FH6)
          FC=FC + FCMXF6 * FMXF6      !          carbon content   (FC6)
          FO=FO + FOMXF6 * FMXF6      !          oxygen content   (FO6)
          FN=FN + FNMXF6 * FMXF6      !          nitrogen content (FN6)
          FNOTAIR=FNOTAIR + FMXF6
        ENDIF
        IF(L0MXF7.NE.0) THEN
          FMXF7=F(L0MXF7+I)        ! mass fraction of material from inlet 7
          FH=FH + FHMXF7 * FMXF7      ! stream * hydrogen content (FH7)
          FC=FC + FCMXF7 * FMXF7      !          carbon content   (FC7)
          FO=FO + FOMXF7 * FMXF7      !          oxygen content   (FO7)
          FN=FN + FNMXF7 * FMXF7      !          nitrogen content (FN7)
          FNOTAIR=FNOTAIR + FMXF7
        ENDIF
        IF(L0FO.NE.0) F(L0FO+I)=FO
        IF(L0FC.NE.0) F(L0FC+I)=FC
        IF(L0FH.NE.0) F(L0FH+I)=FH
        IF(L0FN.NE.0) F(L0FN+I)=FN
        IF(WOOD) THEN                ! -------------------'wood' section
C.... Mass fraction of wood derivatives, FWD:
          FWD=F(L0FWD+I)            ! FWD is treated as a source-free variable
          FH=FH + HINWD*FWD        ! its contribution to h-content of 'gas'
          FC=FC + CINWD*FWD        !              and to c
          FO=FO + OINWD*FWD        !              and to o
          FNOTAIR=FNOTAIR + FWD
          F(L0CHAR+I)=AMIN1(F(L0CHAR+I), FWD*CHINWD) ! char concentration is
          YCHAR=F(L0CHAR+I)                          ! computed elsewhere
          YWOOD=F(L0WOOD+I)                          ! wood also
        ELSE                         ! when 'wood is not present
C.... Coal- (or oil-) burning in gaseous phase (i.e. FGAS=1):
          FGAS=1.0
          YWOOD=0.0
          YCHAR=0.0
          FWD=0.0
        ENDIF                    ! end of 'wood'section
C.... Compute the species concentrations (N2 can be in coal & air):
        FAI=1.0 - FNOTAIR
        FO=FO + 0.232 * FAI
        FN=FN + 0.768 * FAI
        YN2=FN                    ! initial guess
        IF(WOOD) THEN
C.... Mass fraction of phase 1 in the gaseous state:
C     Note that it is necessary to recognise that wood and char may be
C     found anywhere because their burning is kinetically controlled;
          FGAS=1.- YWOOD - YCHAR   ! proportion of phase 1 which is true gas
          RFGAS=1./FGAS             ! its reciprocal
C.... Mass fractions of elements in the gas
          GO=(FO - OINWD*YWOOD)*RFGAS
          GC=(FC - CINWD*YWOOD - YCHAR)*RFGAS
          GH=(FH - HINWD*YWOOD)*RFGAS
        ELSE
          GO =FO
          GC =FC
          GH =FH
        ENDIF
C.... Determine composition of gas on the presumption that mixed is
C     burned:
        ZCO2=3.666667*GC           ! 3.66667=44 / 12=mlwt CO2 / mlwt C
        ZH2O=9.0*GH                ! 9.0   =18 / 2=mlwt H2O / mlwt H
        ZO2=GO - 2.666667*GC - 8.0*GH  ! 2.6667=32 / 12; 8.0=16 / 2
        YVOL=0.0                   !  OK ?????
        IF(ZO2.GE.0.0) THEN         ! region 1: O2, H2O, CO2, N2are present
          ZH2=0.0
          YH2=0.0
          ZCO=0.0
          YCO=0.0
          ZVOL=0.0
          YO2=ZO2*FGAS
          YCO2=ZCO2*FGAS
          YH2O=ZH2O*FGAS
          IF(WOOD) THEN
            YVOL=YVOL + VLINWD*YWOOD   ! ??? how can volatiles be finite on R 1
            YH2O=YH2O + MOINWD*YWOOD
          ENDIF
          YN2=1.0 - YO2 - YCO2 - YH2O   ! ?????
          yn2=amax1(yn2,0.0)
          RMIX=GASCON*(YO2*0.03125 + YH2O*0.055556 + YCO2*0.022727
     1                              + YN2*0.035714 + TINY)
          HSUB=0.0
        ELSE
C.... O2 is not present, i.e. it is region 2 or 3
          ZO2=0.0
          YO2=0.0
C.... Suppose free CO, CO2, H2, H2O exist at the full-oxidation limit,
C     i.e. GOfull=((16/2)*GH + (32/12)*GC)*(1. - GOfull)/(1. - GO)
          GOFULL=(8.*GH + 2.666667*GC)/(1.- GO + 8.*GH + 2.666667*GC)
          FACTFU=(1.- GOFULL)/(1.- GO)
          ZCO2FU=FACTFU*3.666667*GC
          ZH2OFU=FACTFU*9.0*GH
C.... At the part-oxidation limit only CO and H2 exist:
          GOPART=1.333333*GC/(1.- GO + 1.333333*GC)
          FACTPA=(1.- GOPART)/(1.- GO)
          ZCOPA=FACTPA*2.333333*GC
          ZH2PA=FACTPA*GH
C.... For the actual mixture:
          FRAC=(GO - GOPART)/(GOFULL - GOPART)
          IF(FRAC.GE.0.0) THEN
C.... The composition of region 2 (i.e. CO2,H2O,CO,H2 and N2):
            ZVOL=0.0
            YVOL=0.0
            ZCO2=ZCO2FU*FRAC
            ZH2O=ZH2OFU*FRAC
            ZCO=ZCOPA *(1.- FRAC)
            ZH2=ZH2PA *(1.- FRAC)
            YCO2=ZCO2*FGAS
            YH2O=ZH2O*FGAS
            YCO=ZCO*FGAS
            YH2=ZH2*FGAS
            IF(WOOD) THEN
              YVOL=YVOL + VLINWD*YWOOD
              YH2O=YH2O + MOINWD*YWOOD
            ENDIF
            YN2=1.0 - YCO2 - YCO - YH2O - YH2   ! ?????
            yn2=amax1(yn2,0.0)
            RMIX=GASCON*(YCO2*0.022727 + YH2O*0.055556 + YCO*0.035714
     1                                + YH2*0.5 + YN2*0.035714 + TINY)
            HSUB=YCO*HCOCO2 + YH2*HHH2O
          ELSE
C.... The composition of region 3 (i.e. CO,H2,N2 and volatiles):
            ZCO2=0.0
            YCO2=0.0
            ZH2O=0.0
            YH2O=0.0
            ZCO=1.75*GO
            IF(WOOD) THEN
              CONSTA=CINVL
              CONSTB=HINVL
            ELSE
              CONSTA=RHOAG
              CONSTB=RHOBG
            ENDIF
            ZVOL=(GC - ZCO*0.428571)/(CONSTA+tiny)
            ZH2=GH - ZVOL*CONSTB
            YCO=ZCO*FGAS
            YH2=ZH2*FGAS
            YVOL=ZVOL*FGAS
            IF(WOOD) THEN
              YVOL=YVOL + VLINWD*YWOOD
              YH2O=YH2O + MOINWD*YWOOD
            ENDIF
            YN2=1.0 - YCO - YH2O - YH2 - YVOL  ! ?????
            yn2=amax1(yn2,0.0)
            RMIX=GASCON*(YCO*0.035714 + YH2*0.5 + YN2*0.035714
     1                                 + YVOL*RMWTVL + TINY)
            HSUB=YCO*HCOCO2 + YH2*HHH2O + YVOL*HVOL
          ENDIF
        ENDIF
C.... Calculate the absolute temperatures and the first-phase density:
        IF(WOOD) HSUB=HSUB + YCHAR*HCHAR + YWOOD*HWOOD
C... Save HSUB into H0 store
        F(L0H012+I)=HSUB
        TEM=(F(L0H1+I) - HSUB)/F(L0CP1+I)
        TEM=AMIN1(VARMAX(TEMP1), AMAX1(100., TEM, VARMIN(TEMP1)))
        DENS=ABSPREs(I)/(RMIX*TEM)
C.... Wood and char are supposed to contribute their mass to the mixture,
C     without significantly increasing its volume.
        IF(WOOD) DENS=DENS + YWOOD + YCHAR
        DENSITY=AMIN1(VARMAX(DEN1), AMAX1(0., DENS, VARMIN(DEN1)))
        IF(.NOT.ONEPHS) THEN
          TEM2=(F(L0H2+I) - HCHX)/F(L0CP2+I)
C... save HCHX into H02 store
          F(KH02+I)=HCHX
          TEM2=AMIN1(VARMAX(TEMP2), AMAX1(100., TEM2, VARMIN(TEMP2)))
          IF(L0TMP2.NE.0) F(L0TMP2+I)=TEM2
C.... Estimated interface temperature:
          CINT1=CINT(H1)
          CINT2=CINT(H2)
          FACTEM=(TEM*CINT1 + TEM2*CINT2)/(CINT1 + CINT2)
          PHINT1=F(L0CP1+I)*FACTEM + HSUB
          PHINT2=F(L0CP2+I)*FACTEM + HCHX
C.... Interface values used for conductive transport between phases
          F(L0H1IF+I)=PHINT1
          F(L0H2IF+I)=PHINT2
        ENDIF
        IF(LBYCO. NE.0) F(L0YCO+I)=YCO
        IF(LBYCO2.NE.0) F(L0YCO2+I)=YCO2
        IF(LBYH2. NE.0) F(L0YH2+I)=YH2
        IF(LBYH2O.NE.0) F(L0YH2O+I)=YH2O
        IF(LBYO2. NE.0) F(L0YO2+I)=YO2
        IF(LBRMIX.NE.0) F(L0RMIX+I)=RMIX
        IF(LBYN2. NE.0) F(L0YN2+I)=YN2
        IF(LBYVOL.NE.0) F(L0YVOL+I)=YVOL
        IF(l0tem. NE.0) F(L0TEM+I)=TEM
        IF(l0tmp2. NE.0) F(L0TMP2+I)=TEM2
c
      IF(L0HSUB.NE.0) F(L0HSUB+I)=HSUB
c-------------------------------------------- end of 7-gases option
      ELSEIF(IGRND.EQ.8) THEN     ! Density for cvd special-purpose program
        CALL CVDPRP(DENSITY,IFILEP,I,2)
c
      ELSEIF(IGRND.EQ.9) THEN     ! Density from CHEMKIN-database
        CALL GXCHKI(I,DENSITY)
c
      ELSEIF(IGRND.EQ.10) THEN    ! Density=1/linear function of temperature
        IF(L0TEM.GT.0) THEN
          DENSITY=1.0/(RHOAG + RHOBG*(F(L0TEM+I)+TEMP0))
        ELSE
          call writ40('temperature must be solved for or stored')
          call writ1i('igrnd   ',igrnd)
         CALL SET_ERR(215,'Error. Temperature must be solved or stored.'
     *               ,1)
C          call earout(1)
        endif
      ENDIF
      IF(IGRND.NE.-1) THEN
        LBDEN=ITWO(DEN1,DEN2,IPHASE.EQ.1)
        IF(LBDEN.GT.0) THEN
          DENSITY=AMAX1(VARMIN(LBDEN),AMIN1(VARMAX(LBDEN),DENSITY))
          IF(DENSITY.LE.0.0) THEN
            call writ1r('density ',density)
            call writ2i('den1    ',den1,'den2    ',den2)
            call writ1r('density ',density)
            call writ1r('varmin  ',varmin(lbden))
            call writ1r('varmax  ',varmax(lbden))
            call writ40('density.le.0.0, in gxdens; run aborted')
            call writ3i('igrnd   ',igrnd,'cell no.',i,'izstep  ',izstep)
            call writ3i('isweep  ',isweep,'istep   ',istep,
     1                'indvar  ',indvar)
            CALL SET_ERR(216,'Error. See result file.',1)
C            CALL EAROUT(1)
          ENDIF
        ENDIF
      ENDIF
      END
c-----------------------------------------------------------------------
C
      SUBROUTINE SLBDEN(IPILOPT,USEPIL,dbgloc)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/prpcmn'
      COMMON/CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
      COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(21),IPRL,
     1           IBTAU,IGF4(16),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,
     1           IPRPS,IRADX,IRADY,IRADZ,IVFOL
     1      /F0/ IF01(29),L0XYDX,L0XYDY,IF02(3),L0XYRV,L0XYXG,IF03,
     1           L0XYYG,IF04,L0XYDXG,L0XYDYG,IF05(106),L0AHZ,IF06(17),
     1           L0XYDZ,L0XYDZG,IF07(137)
      COMMON/NAMFN/NAMFUN,NAMSUB
cccc      common/denonly/igrnden
      LOGICAL USEPIL,DBGLOC,SLD
      common/dbe/dbear
      logical dbear
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB='SLBDEN'
      if(dbgloc) call banner(1,namsub,220306)
      IGR=9
      ISC=IPROP
      if(dbgloc) call writ1i('ipilopt ',ipilopt)
C.... Call GROUND for the user-set property:
      IF(IPILOPT.EQ.0) THEN
        IF(USEGRD) THEN
          if(dbgloc) then
            call writ40('GROUND is called to set a property...   ')
            call writ2i('Group= ',igr,'Section=',isc)
          endif
          CALL GROUND
        ENDIF
        GO TO 800
C.... For a constant property go to the slab loop
      ELSEIF(IPILOPT.EQ.-1.AND.USEPIL) THEN
        if(dbgloc) call writ40('property is uniform')
        GO TO 700
      ENDIF
c
C.... Set constants and other auxiliary variables:
C-----------------------------------------------------------------------
C.... Density of the first or second phase:
      IF(IPHASE.EQ.1) THEN
        RHOAG=RHO1A
        RHOBG=RHO1B
        RHOCG=RHO1C
      ELSE
        RHOAG=RHO2A
        RHOBG=RHO2B
        RHOCG=RHO2C
      ENDIF
C.... Density using Extended Simple Chemically-Reacting System:
      IF(IPILOPT.EQ.9 .AND. GRNDNO(9,RHOAG)) THEN
        CALL GXSCRI
        GO TO 800
      ENDIF
c.... make necessary initial settings
      CALL SETDEN
c
C.... Loop over slab to get and set cell properties:
  700 IF(IPRPS.EQ.0) THEN ! PRPS is not stored; one material only, no blockages
        IGRND=IPILOPT
        if(dbrho) then
          write(14,*) 'setting phase-',iphase,
     1                ' density via PIL with igrnd=',igrnd
        endif
        DO 60 I=1,NXNYST
  60    F(KPROP+I)=DENSITY(I)   ! set slab values via PIL values
      ELSEIF(USEPIL) THEN
c.... One material only, with blockages
        IGRND=IPILOPT
        DO 70 I=1,NXNYST
          IF(SLD(I)) THEN
            F(KPROP+I)=TINY
          ELSE
            F(KPROP+I)=DENSITY(I)
          ENDIF
  70    CONTINUE
      ELSEIF(SURF .OR. HOL) THEN
C.... Properties are set for either scalar equation method (SURF=T), or
C     height-of-liquid method (HOL=T):
        CALL SURHOL(KPROP,IPROP,IFILEP,IPHASE.EQ.1)
      ELSE
        if(dbgloc) call writ40('prfile called; use imat for each cell')
C.... Properties are set using material flag for an individual cell.
        CALL PRFILE(IPILOPT,dbgloc)
      ENDIF
C----------------------------------------------------------------------
C.... Corrections, debug print-out, and other property adjustments:
        if(dbgloc .and. IPILOPT.eq.7) then
          if(izstep.ge.izdb1 .and.izstep.le.izdb2.and.
     1       isweep.ge.iswdb1.and.isweep.le.iswdb2    ) then
            if(l0mdot.ne.0) call prn('mdt*',-l0mdot)
            if(l0fo.ne.0) call prn('fo  ',-l0fo)
            if(l0fc.ne.0) call prn('fc  ',-l0fc)
            if(l0fh.ne.0) call prn('fh  ',-l0fh)
            if(l0hsub.ne.0) call prn('hsub',-l0hsub)
            if(wood) then
              call prn('wood',-l0wood)
              call prn('char',-l0char)
              call prn('vola',-l0vola)
            endif
            if(lbyco. ne.0) call prn('yco ',lbyco)
            if(lbyco2.ne.0) call prn('yco2',lbyco2)
            if(lbyh2. ne.0) call prn('yh2 ',lbyh2)
            if(lbyh2o.ne.0) call prn('yh2o',lbyh2o)
            if(lbyo2. ne.0) call prn('yo2 ',lbyo2)
            if(lbrmix.ne.0) call prn('rmix',lbrmix)
            if(lbyn2. ne.0) call prn('yn2 ',lbyn2)
            if(lbyvol.ne.0) call prn('yvol',lbyvol)
            if(temp1. ne.0) call prn('tem ',temp1)
            call prn('h1  ',-l0h1)
            if(.not.onephs) then
              call prn('tem2',temp2)
              call prn('h1if',-l0h1if)
              call prn('h2  ',-l0h2)
              call prn('h2if',-l0h2if)
            endif
          endif
        ENDIF
C.... Call GREX to correct a property set above:
 800  IF(USEGRX) CALL GREX3
      IF(USEALT) CALL ALTPRP
      IF(USEGRD) THEN
        IF(IPILOPT.GT.0) CALL GROUND
      ENDIF
      if(dbgloc) call banner(2,'slbden',0)
      NAMSUB='slbden'
      END
c-----------------------------------------------------------------------
C
      SUBROUTINE SETDEN
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/prpcmn'
      COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,NWHOLE,IGSH,
     1   IGF3(37),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
      COMMON /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
      LOGICAL QEQ,LTZ
c-----------------------------------------------------------------------
      IF(SOLVE(1))  L0P=L0F(1)
c???? can the next be simplified?
      IF(IPHASE.EQ.1) THEN
        IF(ITEM1.NE.0) THEN
          L0TEM=L0F(ITEM1)
        ELSEIF(TEMP1.GT.0) THEN
          L0TEM=L0F(TEMP1)
        ELSEIF(TMP1.LE.GRND.OR.QEQ(TMP1,1.11111111)) THEN
          L0TEM=L0F(TEMP1)
        ENDIF
        L0DEN=L0F(DEN1)
      ELSE
        IF(ITEM2.NE.0) THEN
          L0TEM=L0F(ITEM2)
        ELSEIF(TMP2.LE.GRND) THEN
          L0TEM=L0F(TEMP2)
        ENDIF
        L0DEN=L0F(DEN2)
        IF(GEN2.NE.0) L0GEN=L0F(GEN2)
      ENDIF
      L0SCAL=0 ; L0CA=0 ; L0CB=0
      IF(IGRND.EQ.1 .OR. IGRND.EQ.2) THEN
C.... Density is a function of either one variable (enthalpy or other
C     scalar); or two, or three scalars:
        IF(IBUOYB.GT.0) THEN
          IPH=IPHASE - 1
          L0SCAL=L0F(IBUOYB+IPH)
          IF(IBUOYC.GT.0) L0CA=L0F(IBUOYC+IPH)
          IF(IBUOYA.GT.0) THEN
            L0CB=L0F(IBUOYA+IPH)
            RHODG=RHO2A
          ELSE
            L0CB=L0CA
            RHODG=0.0
          ENDIF
        ELSE
          LBH =ITWO(H1,H2,IPHASE.EQ.1)
          LBSCAL=ITWO(LBH , NINT(RHOCG) , QEQ(RHOCG,0.0))
          L0SCAL=L0F(LBSCAL)
        ENDIF
      ELSEIF(IGRND.EQ.3) THEN
C.... Density for the compressible flow and also shallow-water theory:
        SHALWT=LTZ(RHO2) .AND. ONEPHS .AND. IZSTEP.GT.1
        IF(STORE(C3)) L0SCAL=L0F(C3)
      ELSEIF(IGRND.EQ.5) THEN
        IF(STORE(C1)) L0SCAL=L0F(C1 + ITWO(0,1,IPHASE.EQ.1))
      ELSEIF(IGRND.EQ.6) THEN
C.... Ideal-gas law for mixture of three chemical 'species':
        L0FUEL=L0F(LBFUEL) ; L0OXID=L0F(LBOXID) ; L0PROD=L0F(LBPROD)
      ELSEIF(IGRND.EQ.7) THEN
C.... Ideal-gas law for a mixture of the 7 gases involved in a simple
C     model of the combustion in air of one or more hydrocarbons
C     containing C, H, N and an inert ash (coal, oil or wood). The gases
C     are: HCx, O2, N2, CO2, CO, H2O, and H2.
c
c    allowance is made for material entering from up to 8 supply points,
c    the local mass fractions of which being measures by the variables
c    MIXF, MXF1, etc up to MXF7
c
c   Each stream can have its own elemental mass fraction FC, FH, FO, FN
c   which must be read from the Q1 file (except that the values for the
c   MIXF stream are supplied by way of  ????)
        IF(LBMIXF.GT.0) L0MIXF=L0F(LBMIXF)
c
        IF(LBMXF1.GT.0) L0MXF1=L0F(LBMXF1)
        IF(LBMXF2.GT.0) L0MXF2=L0F(LBMXF2)
        IF(LBMXF3.GT.0) L0MXF3=L0F(LBMXF3)
        IF(LBMXF4.GT.0) L0MXF4=L0F(LBMXF4)
        IF(LBMXF5.GT.0) L0MXF5=L0F(LBMXF5)
        IF(LBMXF6.GT.0) L0MXF6=L0F(LBMXF6)
        IF(LBMXF7.GT.0) L0MXF7=L0F(LBMXF7)
c
        IF(LBFC.GT.0) L0FC=    L0F(LBFC)
        IF(LBFH.GT.0) L0FH=    L0F(LBFH)
        IF(LBFO.GT.0) L0FO=    L0F(LBFO)
        IF(LBFN.GT.0) L0FN=    L0F(LBFN)
c
        L0H1 =L0F(H1)
        L0CP1=L0F(ISPH1)
        IF(TEMP1.NE.0) L0TEM=L0F(TEMP1)
        IF(.NOT.ONEPHS) THEN
          L0MDOT=L0F(INTMDT)
          L0H2=L0F(H2)
          L0CP2=L0F(ISPH2)
          IF(TEMP2.NE.0) L0TMP2=L0F(TEMP2)
        ENDIF
        IF(LBYCO.NE.0)  L0YCO=L0F(LBYCO)
        IF(LBYCO2.NE.0) L0YCO2=L0F(LBYCO2)
        IF(LBYH2.NE.0)  L0YH2=L0F(LBYH2)
        IF(LBYH2O.NE.0) L0YH2O=L0F(LBYH2O)
        IF(LBYO2.NE.0)  L0YO2=L0F(LBYO2)
        IF(LBYN2.NE.0)  L0YN2=L0F(LBYN2)
        IF(LBRMIX.NE.0) L0RMIX=L0F(LBRMIX)
        IF(LBYVOL.NE.0) L0YVOL=L0F(LBYVOL)
        IF(WOOD) THEN
          L0WOOD=L0F(LBWOOD)
c moved to start          LBFWD=LBNAME('FWD ')
          L0FWD=L0F(LBFWD)
          L0CHAR=L0F(LBCHAR)
        ENDIF
      ENDIF
cccc      if(flag) call banner(2,'setden',220306)
      END
c-----------------------------------------------------------------------
C
c.... INIDEN is called from INIPRP at the beginning of a run to make
c     once for all settings related to density
c
      SUBROUTINE INIDEN
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/prpcmn'
      COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(39),
     1            ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
     1      /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
      LOGICAL GRNDNO,EQZ,GRPROP(9),GRN,STORE
      DIMENSION NAMPRP(30)
      COMMON/HBASE/IH01,IH02,KH01,KH01H,KH01L,KH02,KH02H,KH02L,L0H012
      EQUIVALENCE (NAMPRP(1),DENST1)
C
      if(dbrho) call banner(1,'INIDEN',220306)
C.... Check what options will be used for first-phase density, and store
c     the result in the grprop array:
      CALL PRPGRN(RHO1,1,GRPROP) ! fill grprop array. Only item 7 is used below
      IF(.NOT.STORE(1)) THEN
        IF(GRNDNO(3,RHO1).OR.GRNDNO(5,RHO1).OR.GRNDNO(6,RHO1).OR.
     1     GRNDNO(7,RHO1).OR.GRNDNO(3,RHO2).OR.GRNDNO(5,RHO2).OR.
     1     GRNDNO(6,RHO2).OR.GRNDNO(7,RHO2) ) THEN
          CALL WRIT40('RHOx=GRND3,5,6,7 requires pressure...  ')
         CALL SET_ERR(217,'Error. RHOx=GRND3,5,6,7 requires pressure...'
     *               ,2)
C          CALL EAROUT(2)
        endif
      endif
C.... Default settings for compressible flow.
      IF(GRNDNO(3,RHO1)) THEN
        DRH1DP=RHO1
      ELSEIF(GRNDNO(3,RHO2)) THEN
        DRH2DP=RHO2
      ELSEIF(GRNDNO(4,RHO1).OR.GRNDNO(5,RHO1).OR.GRNDNO(6,RHO1)) THEN
        IF(ITEM1.EQ.0 .AND. TEMP1.LT.0 .AND. .NOT.GRN(TMP1)) THEN
          CALL WRIT40('RHO1=GRND4,5,6 requires temperature...  ')
        CALL SET_ERR(218,'Error. RHO1=GRND4,5,6 requires temperature...'
     *              ,2)
C          CALL EAROUT(2)
        ENDIF
        IF(GRNDNO(5,RHO1)) DRH1DP=RHO1
      ELSEIF(GRNDNO(4,RHO2).OR.GRNDNO(5,RHO2).OR.GRNDNO(6,RHO2)) THEN
        IF(ITEM2.EQ.0 .AND. TEMP2.LT.0 .AND. .NOT.GRN(TMP2)) THEN
          CALL WRYT40('RHO2=GRND4,5,6 requires temperature...  ')
        CALL SET_ERR(219,'Error. RHO2=GRND4,5,6 requires temperature...'
     *              ,2)
        ENDIF
        IF(GRNDNO(5,RHO2)) DRH2DP=RHO2
      ELSEIF(GRPROP(7) .OR. GRNDNO(7,RHO2)) THEN   !----------- 7 gases model
cccc        LBYCO=LBNAME('YCO ') ; LBYCO2=LBNAME('YCO2')
cccc        LBYH2=LBNAME('YH2 ') ; LBYH2O=LBNAME('YH2O')
cccc        LBYO2=LBNAME('YO2 ')
cccc        LBYN2=LBNAME('YN2 ') ; LBRMIX=LBNAME('RMIX')
cccc        LBYVOL=LBNAME('YVOL')
cccc        LBMIXF1=LBNAME('MXF1') ; LBMXF2=LBNAME('MXF2') ! LBs for
cccc        LBMXF3=LBNAME('MXF3') ; LBMXF4=LBNAME('MXF4') ! up to 7
cccc        LBFC=LBNAME('FC') ; LBFH=LBNAME('FH') ! LBs for
cccc        LBFO=LBNAME('FO') ; LBFN=LBNAME('FN') ! elements C,H,O,N
cccc        LBMXF5=LBNAME('MXF5') ; LBMXF6=LBNAME('MXF6')  ! extra
cccc        LBMXF7=LBNAME('MXF7')                         ! inflows
        IF(LBMXF1.NE.0) THEN
          CALL RQ1R('7GAS','FH1',FHMXF1)
          CALL RQ1R('7GAS','FC1',FCMXF1)
          CALL RQ1R('7GAS','FO1',FOMXF1)
          CALL RQ1R('7GAS','FN1',FNMXF1)
          CALL RQ1R('7GAS','HF1',HMXF1)
        ENDIF
        IF(LBMXF2.NE.0) THEN
          CALL RQ1R('7GAS','FH2',FHMXF2)
          CALL RQ1R('7GAS','FC2',FCMXF2)
          CALL RQ1R('7GAS','FO2',FOMXF2)
          CALL RQ1R('7GAS','FN2',FNMXF2)
          CALL RQ1R('7GAS','HF1',HMXF2)
        ENDIF
        IF(LBMXF3.NE.0) THEN
          CALL RQ1R('7GAS','FH3',FHMXF3)
          CALL RQ1R('7GAS','FC3',FCMXF3)
          CALL RQ1R('7GAS','FO3',FOMXF3)
          CALL RQ1R('7GAS','FN3',FNMXF3)
          CALL RQ1R('7GAS','HF1',HMXF3)
        ENDIF
        IF(LBMXF4.NE.0) THEN
          CALL RQ1R('7GAS','FH4',FHMXF4)
          CALL RQ1R('7GAS','FC4',FCMXF4)
          CALL RQ1R('7GAS','FO4',FOMXF4)
          CALL RQ1R('7GAS','FN4',FNMXF4)
          CALL RQ1R('7GAS','HF1',HMXF4)
        ENDIF
        IF(LBMXF5.NE.0) THEN
          CALL RQ1R('7GAS','FH5',FHMXF5)
          CALL RQ1R('7GAS','FC5',FCMXF5)
          CALL RQ1R('7GAS','FO5',FOMXF5)
          CALL RQ1R('7GAS','FN5',FNMXF5)
          call writ2r('fh5     ',FHMXF5,'fcf     ',fcmxf5)
          call writ2r('fo5     ',FoMXF5,'fnf     ',fnmxf5)
          CALL RQ1R('7GAS','HF1',HMXF5)
        ENDIF
        IF(LBMXF6.NE.0) THEN
          CALL RQ1R('7GAS','FH6',FHMXF6)
          CALL RQ1R('7GAS','FC6',FCMXF6)
          CALL RQ1R('7GAS','FO6',FOMXF6)
          CALL RQ1R('7GAS','FN6',FNMXF6)
          CALL RQ1R('7GAS','HF1',HMXF6)
        ENDIF
        IF(LBMXF7.NE.0) THEN
          CALL RQ1R('7GAS','FH7',FHMXF7)
          CALL RQ1R('7GAS','FC7',FCMXF7)
          CALL RQ1R('7GAS','FO7',FOMXF7)
          CALL RQ1R('7GAS','FN7',FNMXF7)
          CALL RQ1R('7GAS','HF1',HMXF7)
        ENDIF
        IF(GRNDNO(7,RHO1)) THEN
          RHOAG=RHO1A ; RHOBG=RHO1B ; RHOCG=RHO1C
        ELSE
          RHOAG=RHO2A ; RHOBG=RHO2B ; RHOCG=RHO2C
        ENDIF
        IF(.NOT.NULLPR) THEN
          CALL WRITBL
          CALL WRITST
          CALL WRIT40('Initial print-out for RHO1=7GASES model:')
          CALL WRIT3R('CINCL   ',RHOAG,'HINCL   ',RHOBG,
     1                'NINCL   ',RHOCG)
        ENDIF
C.... The heat of combustion per unit mass of CO:
        HCCO2=32.792E6 ; HCCO=9.208E6 ; HHH2O=120.9E6
        HCOCO2=(12.0/28.0) * (HCCO2 - HCCO)
        HCHX=RHOAG*HCCO2 + RHOBG*HHH2O
        HVOL=HCHX
C.... Assume that the volatile gas is propane:
        AMWTVL=44.
        IF(.NOT.NULLPR)
     1  CALL WRIT3R('HCOCO2  ',HCOCO2,'HHH2O   ',HHH2O,'NCHX    ',HCHX)
        IF(.NOT.ONEPHS) THEN
          CALL GXMAKE(L0H1IF,NXNYST,'H1IF')
          CALL GXMAKE(L0H2IF,NXNYST,'H2IF')
        ENDIF
C.... Memory allocation for the debug print-out:
        if(dbrho) then
          call gxmake(l0hsub,nxnyst,'hsub')
          call gxmake(l0vola,nxnyst,'vola')
        endif
C.... Model for the 'wood' burning, i.e. any solid material which shares
c     the gas-phase velocities:
cccc        LBFWD=LBNAME('FWD ')           ! storage of FWD is the signal
        WOOD=LBFWD.NE.0
        IF(.NOT.NULLPR) CALL WRIT1L('WOOD    ',WOOD)
        IF(WOOD) THEN
cccc          LBWOOD=LBNAME('WOOD')        ! but there may be 2 derivatives:
cccc          LBCHAR=LBNAME('CHAR')        ! 'wood' and 'char'
C.... Get data for the wood from SPEDAT:
C.... MOINWD, the moisture fraction of the wood, is defined as the
C     amount of water in unit mass of MOIST wood:
          CALL GETSDR('WOODBURN','MOINWD ',MOINWD) ! read in moisture in wood
C.... CINWD,  HINWD and    VLINWD denote the mass fractions of
c     carbon, hydrogen and volatiles in DRY wood:
          CALL GETSDR('WOODBURN','VLINWD',VLINWD)  ! read in volatile
          CALL GETSDR('WOODBURN','HINVL ',HINVL)   ! and its hydrogen
          CINVL=1.- HINVL                        ! but its carbon
          HINWD=HINVL*VLINWD*(1.- MOINWD) +
     1            2./18.*MOINWD                    ! deduce H in vol and h2o
          CINWD=(1.- HINWD)*(1.- MOINWD)         ! C in wood
          CHINWD=(1.- VLINWD)*(1.- MOINWD)        ! H in wood
          OINWD=16./18.*MOINWD                   ! assume all O in H2O
C.... Heat of combustion:
          HVOL=HHH2O*HINVL + HCCO2*RHOAG
          HCHAR=HCCO2
          HWOOD=HVOL*VLINWD + HCHAR*(1.- VLINWD)
C.... Molecular weight of volatiles:
          CALL GETSDR('WOODBURN','MLWTVL',AMWTVL)
          IF(.NOT.NULLPR) THEN
            CALL WRIT3R('VLINWD  ',VLINWD,'HINVL   ',HINVL,
     1                  'MLWTVL  ',AMWTVL)
            CALL WRIT3R('CINVL   ',RHOAG,'HINWD   ',HINWD,
     1                  'CINWD   ',CINWD )
            CALL WRIT3R('HWOOD   ',HWOOD, 'HVOL    ',HVOL,
     1                  'HCHAR   ',HCHAR)
            CALL WRIT2R('CHENWD  ',CHINWD,'MOINWD  ',MOINWD)
          ENDIF
C.... Store reciprocal of the volatile gas molecular weight:
          RMWTVL=1./(AMWTVL+1.E-10)
        ENDIF    ! ------------------------------end of data input for 'wood'
        IF(.NOT.NULLPR) THEN
          CALL WRITST
          CALL WRITBL
        ENDIF
C.... Note that there is no need to employ a constant specific heat; for
C     coding which allows it to depend upon temperature and composition
C     could also be introduced into this subroutine; but then the heats
C     of combustion would also need to be modified accordingly.
c
C... get index for H0
        L0H012=L0F(ITWO(IH01,IH02,IPHASE.EQ.1))
      ENDIF
C.... Consistency check for for compressibility options:
      IF(GRNDNO(5,DRH1DP)) THEN
        IF(EQZ(RHO1C)) RHO1C=1./1.4
      ENDIF
      IF(GRNDNO(5,DRH2DP)) THEN
        IF(EQZ(RHO2C)) RHO2C=1./1.4
      ENDIF
      if(dbrho) call banner(2,'INIDEN',220306)
      END
c-----------------------------------------------------------------------
c