c

C name htbxgr.for................................................120210
C This file contains four subroutines used for Fan Setting option
C of the HotBox. First is HTBXGR which calculates velocities and
C pressure differences. FANDAT is used for reading fan types from
C the FANDATA file which resides in the local directory. FNAREA
C calculates area for each fan type. FANVAL sets sources for
C each fan type.
C
C This is version allows up to 50 different fan types per simulation.
C********************************************************************
      SUBROUTINE HTBXGR
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      INCLUDE '/phoenics/d_includ/pltcfile'
      PARAMETER (NIG=20, NRG=100, NCG=10)
      PARAMETER (NPRP=50, NVAR=2, NPTS=100)
C NPRP - number of fan types, 50
C NVAR - number of data sets per fan, 2
C NPTS - number of data points per set, 100
      COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15)
     1      ,ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ
     1      ,IVFOL
      COMMON/PRDTA1/PRDATA(NPRP,NVAR,NPTS)
      COMMON/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
      COMMON/CMOBTP/OBJTYP(20)
      COMMON/SWRFAN/RADU1,RADU2,GX0,GY0,GZ0,SWN,ICFAN,IFANN(NPRP)
     1             ,INJDIR
      COMMON/FANMS1/IFANT,VELARR(NPRP),NAMPCH(NPRP)
      DIMENSION IFNARR(NPRP,2),PAREA(NPRP),COSET(NPRP),FLARR(NPRP)
     1         ,FDENA(NPRP)
      DIMENSION FANDP(NPRP),SYSDP(NPRP),LAREA(NPRP),ICHK(NPRP)
      DIMENSION SPDRP(NPRP,1),VFLRT(NPRP),ISCAL(1),DPMAX(NPRP)
      DIMENSION NSLB(NPRP),GEOMF(NPRP),FDENS(NPRP)
      CHARACTER*16 NAMARR(NPRP),TYPARR(NPRP),NAMPCH
      CHARACTER FORM*10,NAMFAN*16,OBJTYP*14,INJDIR*1
      LOGICAL QEQ,QNE,LEZ,LAREA,ISINZZ
      LOGICAL LMCH,DONE,DONEWR,LSYS,MASKPT
      COMMON /L0PVCM/L0PV1(6),NMVRPT,L0VRPT,L0VRAD,NVRINF
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      LOGICAL ZEROED,READFAN
      SAVE
C... If non-VR fan matching or system-curve, call old HTBXGR
      NAMSUB='HTBXGR'
      NAMFUN='NONE  '
      IF(READQ1) THEN
        CALL WRIT40('Old HOTBOX input no longer accepted   ')
        CALL WRIT40('Please run Q1 through VR-Editor before')
        CALL WRIT40('running EARTH again                   ')
        CALL WRYT40('Old HOTBOX input no longer accepted   ')
        CALL WRYT40('Please run Q1 through VR-Editor before')
        CALL WRYT40('running EARTH again                   ')
        CALL SET_ERR(526, 'Error. See result file.',2)
      ENDIF
C... Circular Fans with swirl velocity
      CALL SWIRFAN
      NAMSUB='HTBXGR'
      IXL=IABS(IXL)
C********************************************************************
C
C.... Group 1. Run title and other preliminaries
      IF(IGR.EQ.1) THEN
c............................ SECTION 1 .............................
        IF(ISC.EQ.1) THEN
          IF(.NOT.NULLPR)
     1       CALL WRYT40('Call special Ground HTBXGR.F of: 120210')
C ... Initialise
          IF(NAMGRD.EQ.'HTBX') OBJTYP(5)='APERTURE'
          ISWP=0
          DO 1 I=1,NPRP
            LAREA(I)=.FALSE.
            ICHK(I)=0
            VELARR(I)=0.0
    1     CONTINUE
          CALL SUB2L(LMCH,.FALSE.,LSYS,.FALSE.)
C ... Read Fan data from EARDAT
          CALL GETSDL('FANS','LMTCH',LMCH)
          CALL GETSDL('FANS','LSYST',LSYS)
          IF(LMCH.OR.LSYS) THEN
            CALL WRIT2L('LMCH    ',LMCH,'LSYS    ',LSYS)
            CALL WRYT40('Fan data being read from EARDAT ')
            CALL GETSDI('FANS','FUPDT',IFUPDT)
            IFANT=0
            NAMFAN=' '
C ... Read EARDAT information about each Fan type
            DO IQ1=1,NPRP
              NAMARR(IQ1)=' '
              TYPARR(IQ1)=' '
              NAMFAN(1:3)='FAN'
              IF(IQ1.LT.10) WRITE(NAMFAN(4:4),'(I1)') IQ1
              IF(IQ1.GE.10) WRITE(NAMFAN(4:5),'(I2)') IQ1
              CALL GETSDI(NAMFAN,'VELDIR',IFNARR(IQ1,1))
              CALL GETSDI(NAMFAN,'FANFNC',IFNARR(IQ1,2))
              CALL GETSDC(NAMFAN,'FANNAM',NAMARR(IQ1))
              CALL GETSDC(NAMFAN,'FANTYP',TYPARR(IQ1))
              LL=LENGZZ(TYPARR(IQ1))
              CALL UCASZZ(TYPARR(IQ1),LL)
              CALL GETSDC(NAMFAN,'PATFAN',NAMPCH(IQ1))
              IF(NAMARR(IQ1).NE.' ') IFANT=IFANT+1
            ENDDO
          ENDIF
          IF(LSYS) THEN
            CALL GETSDR('FANS','FLMIN',FLRMIN)
            CALL GETSDR('FANS','FLMAX',FLRMAX)
            CALL GETSDI('FANS','NPOINT',INPT)
          ENDIF
          IF(LMCH) THEN
            CALL WRIT1I('IFUPDT  ',IFUPDT)
          ENDIF
          IF(LSYS) THEN
            CALL WRIT2R('FLRMIN  ',FLRMIN,'FLRMAX  ',FLRMAX)
            CALL WRIT1I('INPT    ',INPT)
          ENDIF
          LHOTDA=1
          ZEROED=.FALSE.
        ENDIF
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
 
      ELSEIF(IGR.EQ.13) THEN
        IF(.NOT.LMCH.AND..NOT.LSYS.AND.STEADY) RETURN
        IF(ISC.EQ.1) THEN
C... Set COefficient for FAN match
          IF(LMCH) THEN
            DO IC=1,IFANT
              IF(NPATCH.EQ.NAMPCH(IC).OR.NPATCH(4:).EQ.NAMPCH(IC)) THEN
                IF((ISTEP.EQ.1.OR.QEQ(PAREA(IC),TINY)).AND.
     1              ISWEEP.EQ.FSWEEP) THEN
C ...  Initialise Fan patch coefficient (jrh: done in Gp19 for restart)
                  IF(QNE(FIINIT(P1),READFI)) THEN
                    PDRP=0.0
                    FLRT=PRINTR(2,1,1,PDRP)
                    VEL=FLRT/(3600.0*0.01)
                    FLRT=0.0
                    PDRP=PRINTR(1,2,IC,FLRT)
                    COSET(IC)=PDRP/(VEL+1.0E-10)
                  ENDIF
C.... Check if volume or mass-flow
                  GEOMF(IC)=1.0
                  IFUNC=IFNARR(IC,2)
C.... Get number of slabs in current patch.
                  CALL GETPTC(NPATCH,RT,I1,I2,J1,J2,K1,K2,IT1,IT2)
                  NSLB(IC)=K2-K1+1
                  IF(MASKPT(IPAT)) THEN
                    DO IVRP=1,NMVRPT
C... loop over object-related patches, to find current patch.
                      KVRPAT=L0VRPT+(IVRP-1)*NVRINF
                      IIP=NINT(F(KVRPAT+1))
C... check if VR patch is current patch
                      IF(IIP.EQ.IPAT) THEN
C... it is, so check geometry factor
                        KVRPAD=L0VRAD+(IVRP-1)*6
                        IOBGFC=NINT(F(KVRPAD+2))
C... if geometry factor is < 19, patch is velocity
                        IF(IOBGFC.LT.19.AND.IFUNC.NE.2) THEN
C.... save area ratio (Area facets/Area cells or 1/Area cells) for edge
                          GEOMF(IC)=F(KVRPAT+4)
                        ENDIF
                      ENDIF
                    ENDDO
                  ENDIF
                ENDIF
C.... Set COefficient for velocity normal to fan
                CALL FN1(CO, COSET(IC))
                RETURN
              ENDIF
            ENDDO
          ENDIF
          RETURN
        ELSEIF(ISC.EQ.12) THEN
          DO IC=1,IFANT
            IF(NAMPCH(IC).EQ.NPATCH.OR.NAMPCH(IC)(4:).EQ.NPATCH) THEN
              IFUNC=IFNARR(IC,2)
              IVELD=IFNARR(IC,1)
              IF(LMCH) THEN
C... Fixflu for P1 for circular Fans
                IF(ICFAN.EQ.2.AND.INDVAR.EQ.P1) then
                  CALL SUB3(L0XG,L0F(XG2D), L0YG,L0F(YG2D)
     1                      ,L0ZG,L0F(ZGNZ))
                  CALL SUB2(L0DEN1,L0F(DEN1), L0VAL,L0F(VAL))
                  DO IX=IXF,IXL
                    DO IY=IYF,IYL
                      I=IY+NY*(IX-1)
c... dis, distance from circle centre to cell centre
                      IF(INJDIR.EQ.'X') THEN
                        DIS=(F(L0YG+I)-GY0)**2+(F(L0ZG+IZ)-GZ0)**2
                      ELSEIF(INJDIR.EQ.'Y') THEN
                        DIS=(F(L0XG+I)-GX0)**2+(F(L0ZG+IZ)-GZ0)**2
                      ELSE
                        DIS=(F(L0XG+I)-GX0)**2+(F(L0YG+I)-GY0)**2
                      ENDIF
                      DIS=SQRT(DIS)
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=F(L0DEN1+I)*VELARR(IC)*IFUNC
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ENDIF
              ELSEIF(LSYS) THEN
C ... Settings for System curve calculations
C ... Calculate Patch area at the beginning only
                IF(ISTEP.EQ.1.AND.ISWEEP.EQ.FSWEEP.AND.INDVAR.EQ.1.AND.
     1            .NOT.LAREA(IC).AND.IZ.GT.IZZ) THEN
                  CALL FNAREA(INTTYP,IFUNC,IXF,IXL,IYF,IYL,AREA,DENA,VL)
                  SYSAR=SYSAR+AREA
                  LAREA(IC)=.TRUE.
                  IZZ=IZ
                ENDIF
C ...
C ... Note: BFC treated as Cartesian
                CALL SUB2(L0VAL,L0F(VAL),L0DEN1,L0F(DEN1))
                L0P1=L0F(P1)
                IF(INDVAR.EQ.1) THEN
                  SETVAL=SYSVEL*IFUNC
C ... Set flow-direction velocity
                ELSEIF(INDVAR.EQ.3.OR.INDVAR.EQ.5.OR.INDVAR.EQ.7) THEN
                  SETVAL=SYSVEL*IVELD
                ENDIF
                CALL SUB3(L0XG,L0F(XG2D), L0YG,L0F(YG2D)
     1                   ,L0ZG,L0F(ZGNZ))
                DO IX=IXF,IXL
                  DO IY=IYF,IYL
                    ICELL=IY+NY*(IX-1)
                    IF(IFUNC.LT.2) THEN
C ... External
                      DPFAN=AMAX1(ABS(F(L0P1+ICELL)),DPFAN)
                    ENDIF
c... dis, distance from Fan circle centre to cell centre
                    INDV=1
                    IF(INJDIR.EQ.'X') THEN
                      INDV=3
                      DIS=(F(L0YG+ICELL)-GY0)**2+(F(L0ZG+IZ)-GZ0)**2
                    ELSEIF(INJDIR.EQ.'Y') THEN
                      INDV=5
                      DIS=(F(L0XG+ICELL)-GX0)**2+(F(L0ZG+IZ)-GZ0)**2
                    ELSE
                      INDV=7
                      DIS=(F(L0XG+ICELL)-GX0)**2+(F(L0YG+ICELL)-GY0)**2
                    ENDIF
                    DIS=SQRT(DIS)
                    IF(INDVAR.EQ.1) THEN
                      F(L0VAL+ICELL)=SETVAL*F(L0DEN1+ICELL)
                    ELSE
                      IF(ICFAN.NE.2.OR.INDV.EQ.INDVAR)
     1                          F(L0VAL+ICELL)=SETVAL
                    ENDIF
                    IF((ICFAN.EQ.2).AND.
     1                 (DIS.LT.RADU1.OR.DIS.GT.RADU2))
     1                 f(l0val+icell)=0.0
                    IF(ICFAN.EQ.2) VELARR(IC)=SYSVEL
                  ENDDO
                ENDDO
              ENDIF
              RETURN
            ENDIF
          ENDDO
C ... Set TM COVALs for linear increase heat source
C ... retained for compatibility with old version
          IF(NAMGRD.EQ.'HTBX'.AND.NPATCH(1:2).EQ.'TM') THEN
            ILEN=0
            IF(ISINZZ(NPATCH(3:3))) ILEN=1
            IF(ILEN.EQ.1.AND.ISINZZ(NPATCH(4:4))) ILEN=2
            IF(ILEN.EQ.2.AND.ISINZZ(NPATCH(5:5))) ILEN=3
            IF(ILEN.EQ.0) GO TO 999
            FORM='(I1)'
            WRITE(FORM(3:3),'(I1)') ILEN
            READ(NPATCH(3:3+ILEN-1),FORM,ERR=999) IADD
            GVAL=RG(40+IADD)
            IFST=IPATNF(9,IREG)
            ILST=IPATNF(10,IREG)
            RCUT=ISTEP-IFST+0.5
            COEF1=GVAL/(ILST-IFST+1)
            SETHS=RCUT*COEF1
            CALL FN1(VAL,SETHS)
          ENDIF
        ELSEIF(ISC.EQ.17.AND.INDVAR.LE.W2) THEN
C ... Set VALues for Fan-Matching
          IF(LMCH.AND.IC.GT.0) THEN
            CALL FANVAL(INTTYP,IFNARR(IC,2),INDVAR,VELARR(IC)
     +                 ,IFNARR(IC,1),FANDP(IC),LAREA(IC)
     +                 ,IXF,IXL,IYF,IYL,PAREA(IC),ICHK(IC),IREG
     +                 ,FDENA(IC))
          ENDIF
        ENDIF
        RETURN
C***************************************************************
C
 
      ELSEIF(IGR.EQ.19) THEN
C ... ------- Section 2 ----- Start of sweep
        IF(ISC.EQ.2) THEN
          IF(LMCH.AND.ISTEP.EQ.1.AND.ISWEEP.EQ.FSWEEP) THEN
C ... Read Fan information from FANDATA file
            IF(.NOT.READFAN) THEN
              READFAN=.TRUE.
              DO IT=1,IFANT
                NAMFAN=TYPARR(IT)
                CALL FANDAT(NAMFAN,IT)
                FLRT=0.0
                DPMAX(IT)=PRINTR(1,2,IT,FLRT)
              ENDDO
            ENDIF
          ENDIF
C ... System curve
          IF(LSYS) THEN
            IF(ISWEEP.EQ.FSWEEP) THEN
              IF(.NOT.DONE) THEN
                CALL WRIT40
     +                 ('System characteristic option activated  ')
                CALL WRITBL
                SYSAR=0.0
                SYSVEL=0.0
                ITRMX=LSWEEP/(INPT+2)
                ITRCUT=ITRMX
                FLRTFX=AMAX1(FLRMIN,1.0)
                FLRTUP=(FLRMAX-FLRMIN)/(INPT+1)
                CALL OPENFL(53)
                WRITE(LUNIT(53),205)
                DONE=LSYS
                DONEWR=.NOT.DONE
                ISCAL(1)=1
                IV=0
              ENDIF
            ELSE
              IF(.NOT.DONEWR)
     +            WRITE(LUNIT(53),206) ISWEEP,DPSYS,FLRTFX
              DONEWR=.NOT.DONEWR
              IF(ISWEEP.EQ.ITRCUT.OR.TOTRES(1).LT.1) THEN
                IV=IV+1
                VFLRT(IV)=FLRTFX
                SPDRP(IV,1)=DPSYS
                IF(IV.LT.INPT+2) THEN
                  ITRCUT=ITRCUT+ITRMX
                  FLRTFX=FLRTFX+FLRTUP
                  IF(QEQ(FLRTFX,FLRMAX)) LSWEEP=ISWEEP+ITRMX
                ENDIF
              ENDIF
            ENDIF
            IF(SYSAR.GT.0.0) THEN
              SYSVEL=FLRTFX/(3600.0*SYSAR)
            ELSE
              SYSVEL=0.0
            ENDIF
          ELSE
C ... Fan Matching
            IF(LMCH) THEN
C ... Check that for unsteady fan area has been set
              IF(ISTEP.GT.1) THEN
                DO 1923 IT=1,IFANT
                  PAREA(IT)=AMAX1(PAREA(IT),TINY)
 1923           CONTINUE
              ENDIF
              IF(.NOT.DONE) THEN
                CALL WRITST
                CALL WRIT40('Fan Matching option activated           ')
                DONE=LMCH
              ENDIF
              IF((ISWEEP.GT.1.OR.QEQ(FIINIT(P1),READFI)).AND.
     1            MOD(ISWEEP,IFUPDT).EQ.0) THEN
                IF(ISTEP.EQ.FSTEP.AND.ISWEEP.EQ.FSWEEP.AND.
     1             QEQ(FIINIT(P1),READFI)) THEN
C ... Set source to zero in the first sweep of the restart
                  IF(.NOT.ZEROED) THEN
                    DO 1924 IT=1,IFANT
                      FANDP(IT)=0.0
                      COSET(IT)=1E6
                      VELARR(IT)=0.0
 1924               CONTINUE
                    ZEROED=.TRUE.
                  ENDIF
                ELSE
                  DO 192 IT=1,IFANT
                    PDRP=SYSDP(IT)
                    IF(PDRP.GT.DPMAX(IT)) THEN
                      FLARR(IT)=FLARR(IT)/2.0
                    ELSE
                      FLARR(IT)=PRINTR(2,1,IT,PDRP)
                    ENDIF
                    IFUNC=IFNARR(IT,2)
C ... Calculate slope at the point to linearize momentum source
C.... GEOMF =1 for internal fan, = Area facets/Area cells for fan at
C.... domain edge.
                    VELNOM=FLARR(IT)/(3600.0*PAREA(IT)*GEOMF(IT))
                    FLRPL1=FLARR(IT)+1.0
                    VELPL1=FLRPL1/(3600.0*PAREA(IT)*GEOMF(IT))
                    PRPL1=PRINTR(1,2,IT,FLRPL1)
                    FLRM1=FLARR(IT)-1.0
                    VELM1=FLRM1/(3600.0*PAREA(IT)*GEOMF(IT))
                    PRM1=PRINTR(1,2,IT,FLRM1)
                    SLOPE=(PRM1-PRPL1)/(VELPL1-VELM1)
                    IF(IFANN(it).NE.2.or.IFNARR(IT,2).EQ.2) THEN
C... Internal fan, or rectangular fan at domain edge
                      COSET(IT)=SLOPE*20.0
                      VELVAL=VELNOM+PDRP/SLOPE/20.0
                      VELARR(IT)=VELVAL
                    ELSE
                      COSET(IT)=ONLYMS
                      VELARR(IT)=VELARR(IT)-(VELNOM-VELARR(IT))*
     1                                                    (PDRP/SLOPE)
                    ENDIF
  192             CONTINUE
                ENDIF
              ENDIF
  205         FORMAT(1X,4HISWP,9X,3HSDP,9X,4HFLRT)
  206         FORMAT(1X,I4,1X,1PE11.3,2X,1PE11.3)
            ENDIF
          ENDIF
C ... ------- Section 6 ----- Finish of iz slab
        ELSEIF(ISC.EQ.6) THEN
          IF(ISWEEP.LE.FSWEEP+1) THEN
            DO IT=1,NPRP
              LAREA(IT)=.FALSE.
            ENDDO
          ELSE
            DO IT=1,NPRP
              LAREA(IT)=.TRUE.
            ENDDO
          ENDIF
C ... Internal fan needs loop
          IF(LSYS.AND.IFUNC.EQ.2) THEN
            CALL HILOXY(L0F(P1),P1MX,P1MN,ILMX,ILMN,.false.)
            PHIGH=AMAX1(PHIGH,P1MX)
            IF(IZ.EQ.1) PLOW=PHIGH
            PLOW=AMIN1(PLOW,P1MN)
          ENDIF
          IF((ISWEEP.EQ.LSWEEP.OR.ENUFSW).AND.ITEM1.GT.0) THEN
            IF(IZ.EQ.1) THIGH=-GREAT
            CALL HILOXY(L0F(ITEM1),TM1MX,TM1MN,ILMX,ILMN,.FALSE.)
            IF(TM1MX.GE.THIGH) THEN
              IYLCMX=MOD(ILMX-1,NY)+1
              IXLCMX=(ILMX-IYLCMX)/NY+1
              IZLCMX=IZ
              THIGH=TM1MX
            ENDIF
            IF(IZ.EQ.1) TLOW=THIGH
            TLOW=AMIN1(TLOW,TM1MN)
          ENDIF
C   * ------------------- SECTION 7 ---- Finish of sweep.
        ELSEIF(ISC.EQ.7) THEN
          IF(LMCH.OR.LSYS) THEN
C... Pressure drop through the system
            IF(LMCH) THEN
              DO 1974 IT=1,IFANT
                IF(NEZ(FDENA(IT))) THEN
                  SYSDP(IT)=FANDP(IT)
                  FDENS(IT)=FDENA(IT)
                ENDIF
                IF(ISWEEP.NE.LSWEEP.AND..NOT.ENUFSW) THEN
                  FANDP(IT)=0.0
                  FDENA(IT)=0.0
                ENDIF
 1974         CONTINUE
            ELSE
              DPINT=PHIGH-PLOW
              DPFAN=AMAX1(ABS(DPINT),DPFAN)
              PHIGH=PLOW
              DPSYS=DPFAN
              DPFAN=0.0
            ENDIF
            IF(LMCH.AND.(ISWEEP.EQ.LSWEEP.OR.ENUFSW)) THEN
C ... Overall system flow-rate (source of R1)
              FLRTS=0.0
              DO 1971 IR=1,NUMREG
                CALL GETSO(IR,R1,SORCE)
                IF(SORCE.GT.0.0) FLRTS=FLRTS+SORCE
 1971         CONTINUE
              CALL WRITST
              CALL WRITBL
              CALL WRIT40('Final result: Operating point data      ')
              IF(.NOT.STEADY) THEN
                CALL WRIT1I('TIME STP',ISTEP)
              ENDIF
              IF(LEZ(BUOYD)) BUOYD=1.189
              CALL WRITBL
C... FL 6/12/99 No System flow rate if no external fans
              IFUNC1=2
              DO IT=1,IFANT
                IF(IFNARR(IT,2).NE.2) IFUNC1=0
              ENDDO
              IF(FLRTS.GT.0.0) THEN
                WRITE(LUPR1,'(1X,1A,1P1E12.4,1A)')
     1               'System mass flow rate: ',FLRTS*3600.0
     1              ,' kg/h'
              ELSE
                WRITE(LUPR1,'(1X,1A,1A)')
     1               'The following are internal Fans only,'
     1              ,'no system mass flow rate available'
              ENDIF
              CALL WRITBL
              DO 1972 IT=1,IFANT
                WRITE(LUPR1,'(1X,4A)')
     1               'Fan Name: ',NAMARR(IT),' Fan Type: ',TYPARR(IT)
                FANVEL=FLARR(IT)/(3600.0*PAREA(IT)*GEOMF(IT))
                FANFLR=FLARR(IT)*FDENS(IT)/NSLB(IT)
                WRITE(LUPR1,'(1X,2A,1P1E12.4,1A)')
     1             NAMARR(IT),' mass flow rate:   ',FANFLR
     1            ,' kg/h'
                WRITE(LUPR1,'(1X,2A,1P1E12.4,1A)')
     1             NAMARR(IT),' volume flow rate: ',FLARR(IT)
     1            ,' m^3/h'
                WRITE(LUPR1,'(1X,2A,1P1E12.4,1A)')
     1             NAMARR(IT),' velocity:         ',ABS(FANVEL)
     1            ,' m/s'
                WRITE(LUPR1,'(1X,2A,1P1E12.4,1A)')
     1             NAMARR(IT),' pressure drop:    ',SYSDP(IT)
     1            ,' Pa'
                CALL WRITBL
 1972         CONTINUE
              IF(.NOT.STEADY) THEN
                CALL WRITBL
              ENDIF
              CALL WRITST
            ENDIF
          ENDIF
          IF(ISWEEP.EQ.LSWEEP.OR.ENUFSW) THEN
C... print sweep
            IF(ITEM1.GT.0) THEN
              CALL WRITBL
              CALL WRIT40('Maximum and minimum temperatures:       ')
              IF(.NOT.STEADY) CALL WRIT1I('TIME STP',ISTEP)
              CALL WRIT2R('MAX.TEMP',THIGH,'MIN.TEMP',TLOW)
              CALL WRITBL
              CALL WRIT40('Location of the maximum-temperature cell')
              CALL WRIT3I('X-LOC   ',IXLCMX,'Y-LOC   ',IYLCMX,
     1                    'Z-LOC   ',IZLCMX)
C... Put max.temp cell into plotting file 'Hotspot'
              CALL CLOSFL(52)
              NMFIL(52)='HOTSPOT'
              CALL OPENFL(52)
              LUG=LUNIT(52)
C... Write message
              WRITE(LUG,'(1X,A,i3,a,i3,a,i3)')
     1     'MSG Hot spot location: x=',ixlcmx,' y=',iylcmx,' z=',izlcmx
              WRITE(LUG,'(1X,A,1PE12.4,a)')
     1     'MSG    Max. temperature ',THIGH,' deg'
C... Volume
              WRITE(LUG,1912) 'X',IXLCMX+1,'Y',IYLCMX,IYLCMX,'Z',IZLCMX,
     1          IZLCMX
              WRITE(LUG,1912) 'X',IXLCMX,'Y',IYLCMX,IYLCMX,'Z',IZLCMX,
     1          IZLCMX
              WRITE(LUG,1912) 'Y',IYLCMX+1,'X',IXLCMX,IXLCMX,'Z',IZLCMX,
     1          IZLCMX
              WRITE(LUG,1912) 'Y',IYLCMX,'X',IXLCMX,IXLCMX,'Z',IZLCMX,
     1          IZLCMX
 1912         FORMAT(1X,'GR OU ',A,I6,2(1X,A,2I6),' COL 15')
              CALL CLOSFL(52)
            ENDIF
C...
            IF(LSYS) THEN
              CALL CLOSFL(52)
              NMFIL(52)='SYSDAT'
              CALL OPENFL(52)
              LUG=LUNIT(52)
              IF(LHOTDA.EQ.1)WRITE(LUG,'(4X,1A,10X,1A)')'SVOL','SDP'
              LHOTDA=2
              DO IV=1,INPT+2
                IF(VFLRT(IV).GT.0.0.AND.SPDRP(IV,1).GT.0.0) THEN
                  WRITE(LUG,'(1X,1P2E12.3)')VFLRT(IV),SPDRP(IV,1)
                ENDIF
              ENDDO
              CALL CLOSFL(52)
              CALL CLOSFL(53)
              CALL WRITBL
              CALL WRIT40('System characteristic points during the ')
              CALL WRIT40('course of calculation:                  ')
              CALL GRAPH(VFLRT,50,INPT+2,'SVOL',SPDRP,1,1,
     1                  'SDP ',ISCAL,0.5,0.5,3,LUPR1)
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      RETURN
C
C... Error exit
  999 CALL WRIT40('Error in number denoting a TM patch.    ')
      CALL WRIT40('Check Q1 for TMx, where x must be 1 - 99')
      CALL SET_ERR(323, 'Error in number denoting a TM patch.',1)
C      CALL EAROUT(1)
      END
C***********************************************************************
      SUBROUTINE FANDAT(NAMFAN,ITYP)
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/pltcfile'
      PARAMETER (NPRP=50, NVAR=2, NPTS=100)
c nprp - number of fans 50
c nvar - number of data sets per fan 2
c npts - number of data points per set 100
      COMMON/PRDTA1/PRDATA(NPRP,NVAR,NPTS)
      CHARACTER*256 NMSAV
      CHARACTER*40 LINE
      CHARACTER*16 NAMFAN
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
c
      NAMFUN='FANDAT'
      IREF=54
      IERR=-1
      NMSAV=NMFIL(IREF)
      NMFIL(IREF)='fandata'
C First attempt to read user fandata file in local directory
      CALL OPENZZ(IREF,IERR)
      IF(IERR.NE.0) THEN
C If error, attempt to read default fandata file
        NMFIL(IREF)=NMSAV
        CALL OPENZZ(IREF,IERR)
        IF(IERR.NE.0) THEN
          WRITE(LUPR1,'(a)') ' Cannot open fan data file! in '
          WRITE(LUPR1,'(a)') NMFIL(IREF)
          WRITE(LUPR3,'(a)') ' Cannot open fan data file! in '
          WRITE(LUPR3,'(a)') NMFIL(IREF)
          CALL SET_ERR(324, 'Error. Cannot open fan data file.',2)
        ENDIF
      ENDIF
      L1=LENGZZ(NMFIL(IREF))
      L2=LENGZZ(NAMFAN)
      WRITE(LUPR1,'(A,A)') ' fan data read from file ',NMFIL(IREF)(1:L1)
      WRITE(LUPR1,'(A,A)') ' for fan type ',NAMFAN(1:L2)
      WRITE(LUPR3,'(A,A)') ' fan data read from file ',NMFIL(IREF)(1:L1)
      WRITE(LUPR3,'(A,A)') ' for fan type ',NAMFAN(1:L2)
      LUFD=LUNIT(IREF)
   10 READ(LUFD,'(A)',ERR=100,END=100) LINE
C... Upper-case line to avoid problems with fan names
      LL=LENGZZ(LINE)
      CALL UCASZZ(LINE,LL)
      IF(INDEX(LINE,NAMFAN).EQ.0) GOTO 10
      READ(LUFD,*) INUM
      PRDATA(ITYP,1,1)=INUM
      PRDATA(ITYP,2,1)=INUM
      DO 20 I=1,INUM
         READ(LUFD,*) RVL,RDP
         PRDATA(ITYP,1,I+1)=RVL
         PRDATA(ITYP,2,I+1)=RDP
   20 CONTINUE
  100 IF(INDEX(LINE,NAMFAN).EQ.0) GOTO 110
      CALL CLOSFL(IREF)
      RETURN
  110 CALL WRIT40('Could not find specified fan in FANDATA!')
      WRITE(LUPR1,'(1X,A,1X,A16)') 'Requested name was',NAMFAN
      CALL SET_ERR(325,
     *      'Error. Could not find specified fan in FANDATA.',1)
      END
C***********************************************************************
      SUBROUTINE FNAREA(ITYP,IFNC,IXFF,IXLL,IYFF,IYLL,AREA,DENA,VEL)
C ... This subroutine determines Fan area to be used for calculating
C ... mass-flow condition which is read from the fan-curve.
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      PARAMETER (NPRP=50)
      COMMON/SWRFAN/RADU1,RADU2,GX0,GY0,GZ0,SWN,ICFAN,IFANN(NPRP)
     1             ,INJDIR
      CHARACTER INJDIR*1
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      LOGICAL MASKPT,L2DVRO
C ...
      NAMFUN='FNAREA'
C ... Determine type for Area calculation
      IF(ITYP.EQ.2) THEN ! East
        IF(BFC) THEN
          L0AR=L0B(5)
        ELSE
          L0AR=L0F(AEAST)
        ENDIF
        L0VEL=L0F(U1)
        IF(IXFF.EQ.NX) L0VAL=L0F(WEST(U1))
      ELSEIF(ITYP.EQ.3) THEN ! West
        IF(BFC) THEN
          L0AR=L0B(6)
        ELSE
          L0AR=L0F(AEAST)
        ENDIF
        L0VEL=L0F(WEST(U1))
        IF(IXFF.EQ.1) L0VEL=L0F(U1)
      ELSEIF(ITYP.EQ.4) THEN ! North
        IF(BFC) THEN
          L0AR=L0B(7)
        ELSE
          L0AR=L0F(ANORTH)
        ENDIF
        L0VEL=L0F(V1)
        IF(IYFF.EQ.NY) L0VAL=L0F(SOUTH(V1))
      ELSEIF(ITYP.EQ.5) THEN ! South
        IF(BFC) THEN
          L0AR=L0B(8)
        ELSE
          L0AR=L0F(ANORTH)
        ENDIF
        L0VEL=L0F(SOUTH(V1))
        IF(IYFF.EQ.1) L0VEL=L0F(V1)
      ELSEIF(ITYP.EQ.6) THEN ! High
        IF(BFC) THEN
          L0AR=L0B(9)
        ELSE
          L0AR=L0F(AHIGH)
        ENDIF
        L0VEL=L0F(W1)
        IF(IZ.EQ.NZ) L0VEL=L0F(LOW(W1))
      ELSEIF(ITYP.EQ.7) THEN ! High
        IF(BFC) THEN
          L0AR=L0B(10)
        ELSE
          L0AR=L0F(AHIGH)
        ENDIF
        L0VEL=L0F(LOW(W1))
        IF(IZ.EQ.1) L0VEL=L0F(W1)
      ENDIF
      AREA=0.0
      CALL SUB3(L0XG,L0F(XG2D), L0YG,L0F(YG2D), L0ZG,L0F(ZGNZ))
      DENA=0.0
      TVOL=0.0
      IPW=0
      L0DEN1=L0F(DEN1)
      DO 1332 IX=IXFF,IXLL
        DO 1332 IY=IYFF,IYLL
          IF(MASKPT(IPAT)) THEN ! skip cells not covered by object
            IPW=IPW+1
            IF(.NOT.L2DVRO(IPW)) GO TO 1332
          ENDIF
          I=IY+NY*(IX-1)
          IF(ICFAN.EQ.2) THEN
            IF(INJDIR.EQ.'X') THEN
              DIS=(F(L0YG+I)-GY0)**2+(F(L0ZG+IZ)-GZ0)**2
            ELSEIF(INJDIR.EQ.'Y') THEN
              DIS=(F(L0XG+I)-GX0)**2+(F(L0ZG+IZ)-GZ0)**2
            ELSE
              DIS=(F(L0XG+I)-GX0)**2+(F(L0YG+I)-GY0)**2
            ENDIF
            DIS=SQRT(DIS)
            IF(DIS.LT.RADU1.OR.DIS.GT.RADU2) GOTO 1332
          ENDIF
          IF(F(L0DEN1+I).LE.2.0) THEN
            AREA=AREA+F(L0AR+I)             ! sum area
            DENA=DENA+F(L0DEN1+I)*F(L0AR+I) ! sum dens*area
            TVOL=TVOL+F(L0AR+I)*F(L0VEL+I)  ! sum vel*area
          ENDIF
 1332 CONTINUE
      DENA=DENA/(AREA+TINY) ! average area
      VEL=TVOL/(AREA+TINY)  ! average velocity
      END
C***********************************************************************
      SUBROUTINE FANVAL(ITYP,IFNC,IVAR,FVEL,IDIRV,FDP,LDONE,IXFF,IXLL,
     1                  IYFF,IYLL,FAREA,ICHK,IR,PDENA)
C ... This subroutine calculates and sets value for each Fan type
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      PARAMETER (NPRP=50)
      COMMON/SWRFAN/RADU1,RADU2,GX0,GY0,GZ0,SWN,ICFAN,IFANN(NPRP)
     1             ,INJDIR
      CHARACTER INJDIR*1
      LOGICAL LDONE, DONE1, LSETV, EQZ, QEQ
      SAVE LBNMUC, LBNMVC, LBNMWC
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
c
      NAMFUN='FANVAL'
      DONE1 = ISWEEP .GT. FSWEEP
      LSETV = .NOT.BFC.OR.IFNC.EQ.2
      IF(BFC.AND..NOT.DONE1) THEN
C... Find pointers for 1st phase cartesian components
        LBNMUC=LBNAME('UCRT')
        LBNMVC=LBNAME('VCRT')
        LBNMWC=LBNAME('WCRT')
      ENDIF
C ... Calculate Patch area at the beginning only
      IF(ISWEEP.EQ.FSWEEP.AND..NOT.LDONE) THEN
        IF(ICHK.EQ.0.OR.ICHK.EQ.IR) THEN
          CALL FNAREA(ITYP,IFNC,IXFF,IXLL,IYFF,IYLL,AREA,DENA,VEL)
          FAREA=FAREA+AREA
C... use recovered velocity for restart, else start from guessed 0.1
          IF(QEQ(FIINIT(P1),READFI)) THEN
            FVEL=VEL
          ELSE
            FVEL=0.01
          ENDIF
          LDONE=.TRUE.
          ICHK=IR
        ENDIF
      ENDIF
C ... Calculate Product of Patch area and Density
      IF(LDONE) THEN
        IF(ICHK.EQ.0.OR.ICHK.EQ.IR) THEN
          CALL FNAREA(ITYP,IFNC,IXFF,IXLL,IYFF,IYLL,AREA,DENA,VEL)
          PDENA=PDENA+DENA
          LDONE=.FALSE.
          ICHK=IR
        ENDIF
      ENDIF
C ... Fan Setting part
      CALL SUB2(L0P1,L0F(P1),L0VAL,L0F(VAL))
      L0NGBR=L0F(P1)
      IF(ITYP.EQ.7) THEN
        L0P1=L0F(HIGH(P1))
      ELSEIF(ITYP.EQ.3) THEN
        L0P1=L0F(EAST(P1))
      ELSEIF(ITYP.EQ.5) THEN
        L0P1=L0F(NORTH(P1))
      ELSEIF(IFNC.EQ.2) THEN
        IF(ITYP.EQ.2) THEN
          L0NGBR=L0F(EAST(P1))
        ELSEIF(ITYP.EQ.4) THEN
          L0NGBR=L0F(NORTH(P1))
        ELSEIF(ITYP.EQ.6) THEN
          L0NGBR=L0F(HIGH(P1))
        ENDIF
      ENDIF
C ... Set flow-direction velocity
 
      IF(IVAR.GE.3.OR.IVAR.LE.7) SETVAL=ABS(FVEL)*IDIRV
C ... BFC part
      IF(.NOT.LSETV) THEN
        CALL SUB3R(UCRT,0.0,VCRT,0.0,WCRT,0.0)
        IUVW=3*((IVAR-U1)/2)
        IF(IVAR.EQ.3) THEN
          UCRT=SETVAL
        ELSEIF(IVAR.EQ.5) THEN
          VCRT=SETVAL
        ELSEIF(IVAR.EQ.7) THEN
          WCRT=SETVAL
        ENDIF
C... Calculate inlet velocity component over current PATCH
        CALL BFCVFX(IUVW+1,IUVW+2,IUVW+3,
     1              UCRT,VCRT,WCRT,IZSTEP,VAL,EASP1,NX,NY)
      ENDIF
      CALL SUB3(L0XG,L0F(XG2D), L0YG,L0F(YG2D), L0ZG,L0F(ZGNZ))
      DO 1351 IX=IXFF,IXLL
      DO 1351 IY=IYFF,IYLL
        ICELL=IY+NY*(IX-1)
        IPCL=ICELL
C ... Fan pressure drop different for internal and external
        IF(IFNC.EQ.2) THEN
C ... Internal
          DPINT=F(L0P1+ICELL)-F(L0NGBR+ICELL)
          FDP=AMAX1(ABS(DPINT),FDP)
        ELSE
C ... External
          IF(ITYP.EQ.3) IPCL=IY+NY*IX
          IF(ITYP.EQ.5) IPCL=IY+1+NY*(IX-1)
          FDP=AMAX1(ABS(F(L0P1+IPCL)),FDP)
        ENDIF
        IF(ICFAN.EQ.2) THEN
          IF(INJDIR.EQ.'X') THEN
            DIS=(F(L0YG+ICELL)-GY0)**2+(F(L0ZG+IZ)-GZ0)**2
          ELSEIF(INJDIR.EQ.'Y') THEN
            DIS=(F(L0XG+ICELL)-GX0)**2+(F(L0ZG+IZ)-GZ0)**2
          ELSE
            DIS=(F(L0XG+ICELL)-GX0)**2+(F(L0YG+ICELL)-GY0)**2
          ENDIF
          DIS=SQRT(DIS)
          IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
            IF(LSETV) F(L0VAL+ICELL)=SETVAL
          ELSE
            IF(LSETV) F(L0VAL+ICELL)=0.0
          ENDIF
        ELSE
          IF(LSETV) F(L0VAL+ICELL)=SETVAL
        ENDIF
 1351 CONTINUE
      END
c