C

C name htbxgr.for................................................200121
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 'farray'
      INCLUDE 'd_earth/parvar'
      INCLUDE 'patcmn'
      INCLUDE 'parear'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      INCLUDE '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 /CGLDIM/NX0,NY0,NZ0
      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),TURBN(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),VELFAN(NPRP),
     1          FANVEL(NPRP),TMAX(NPROC),ITMAX(NPROC),IPROCM(1)
      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,DONE_SUM
      COMMON /L0PVCM/L0PV1(6),NMVRPT,L0VRPT,L0VRAD,NVRINF
      COMMON /INDAUX/L0DUM(16),L0TTFLM,L0DUM2(3)
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      LOGICAL ZEROED,READFAN,L2DVRO,NEZ,LPAR,ISACTIVECELLS,LOG1,dbg
      REAL RDUM(6)
      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==1) THEN
c............................ SECTION 1 .............................
        IF(ISC==1) THEN
          IF(.NOT.NULLPR) THEN
            CALL WRYT40('Call special Ground HTBXGR.F of: 200121')
            CALL WRIT40('Call special Ground HTBXGR.F of: 200121')
          ENDIF
C ... Initialise
          IF(NAMGRD=='HTBX') OBJTYP(5)='APERTURE'
          ISWP=0
          LPAR=MIMD.AND.NPROC>1
          IFREQ=IDISPA
          IF(IFREQ==0) IFREQ=LSWEEP+1
          DO 1 I=1,NPRP
            LAREA(I)=.FALSE.
            ICHK(I)=0
            VELARR(I)=0.0
            GEOMF(I)=1.0
    1     CONTINUE
          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)
            RELX=0.5; CALL GETSDR('FANS','RELAX',RELX)
            dbg=.false.; call GETSDL('FANS','DEBUG',dbg)
            IFANT=0
            NAMFAN=' '
C ... Read EARDAT information about each Fan type
            DO IQ1=1,NPRP
              NAMARR(IQ1)=' '
              TYPARR(IQ1)=' '
              NAMFAN(1:3)='FAN'
              TURBN(IQ1)=0.0
              IF(IQ1<10) WRITE(NAMFAN(4:4),'(I1)') IQ1
              IF(IQ1>=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)/=' ') IFANT=IFANT+1
              CALL GETSDR(NAMFAN,'TURBN',TURBN(IQ1))
            ENDDO
          ENDIF
          IF(LSYS) THEN
            CALL GETSDR('FANS','FLMIN',FLRMIN)
            CALL GETSDR('FANS','FLMAX',FLRMAX)
            CALL GETSDI('FANS','NPOINT',INPT)
            CALL WRIT2R('FLRMIN  ',FLRMIN,'FLRMAX  ',FLRMAX)
            CALL WRIT1I('INPT    ',INPT)
            DONE_SUM=.FALSE.; ITRCUT0=0
          ENDIF
          IF(LMCH) THEN
            CALL WRIT1I('IFUPDT  ',IFUPDT)
          ENDIF
          LHOTDA=1
          ZEROED=.FALSE.
        ENDIF
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
 
      ELSEIF(IGR==13) THEN
        IF(.NOT.LMCH.AND..NOT.LSYS.AND.STEADY) RETURN
        IF(ISC==1) THEN
C... Set COefficient for FAN match
          IF(LMCH) THEN
            DO IC=1,IFANT
              IF(NPATCH==NAMPCH(IC).OR.NPATCH(4:)==NAMPCH(IC)) THEN
                IF((ISTEP==1.OR.QEQ(PAREA(IC),TINY)).AND.
     1              ISWEEP==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)
                  IF(MASKPT(IPAT)) THEN ! facetted object
                    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==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<19.AND.IFUNC/=2) THEN
C.... save area ratio (Area facets/Area cells or 1/Area cells) for edge
                          GEOMF(IC)=F(KVRPAT+4)
                        ENDIF
                        IPW=0; NCEL=0
                        IF(IZ==K1) NSLB(IC)=0
                        DO IXX=I1,I2; DO IYY=J1,J2 ! loop over cells on slab
                          IPW=IPW+1
                          IF(.NOT.ISACTIVECELLS(IXX,IYY,IZ)) CYCLE
                          IF(L2DVRO(IPW)) NCEL=NCEL+1 ! count cells in object
                        ENDDO; ENDDO
                        IF(NCEL>0) NSLB(IC)=NSLB(IC)+1 ! if any cells, increment
                      ENDIF
                    ENDDO
                  ELSE ! non-faceted object
                    IF(IZ==K1) NSLB(IC)=0
                    NCEL=0
                    DO IXX=I1,I2; DO IYY=J1,J2 ! loop over cells on slab
                      IF(.NOT.ISACTIVECELLS(IXX,IYY,IZ)) CYCLE
                      NCEL=NCEL+1 ! count cells in object
                    ENDDO; ENDDO
                    IF(NCEL>0) NSLB(IC)=NSLB(IC)+1 ! if any cells, increment
                  ENDIF
                ENDIF
C.... Set COefficient for velocity normal to fan
                CALL FN1(CO, COSET(IC))
                RETURN
              ENDIF
            ENDDO
          ENDIF
          RETURN
        ELSEIF(ISC==12) THEN
          DO IC=1,IFANT
            IF(NAMPCH(IC)==NPATCH.OR.NAMPCH(IC)(4:)==NPATCH) THEN
              IFUNC=IFNARR(IC,2)
              IVELD=IFNARR(IC,1)
              IF(LMCH.AND.ICFAN==2) THEN
C... Fixflu for P1 for circular Fans
                IF(INDVAR==P1) THEN
                  L0XG=L0F(XG2D); L0YG=L0F(YG2D); L0ZG=L0F(ZGNZ)
                  L0DEN1=L0F(DEN1); L0VAL=L0F(VAL); IPW=0
                  DO IX=IXF,IXL
                    DO IY=IYF,IYL
                      I=IY+NY*(IX-1)
                      IF(MASKPT(IPAT)) THEN
                        IPW=IPW+1
                        IF(.NOT.L2DVRO(IPW)) THEN
                          F(L0VAL+I)=0.0; CYCLE
                        ENDIF
                      ENDIF
                        F(L0VAL+I)=F(L0DEN1+I)*VELARR(IC)*IFUNC
                    ENDDO
                  ENDDO
                ELSEIF(INDVAR==KE.OR.INDVAR==EP) THEN
                  ENINL=MAX(TURBN(IC)*VELARR(IC),0.001)
                  IF(INDVAR==KE) THEN
                    CALL FN1(VAL,ENINL)
                  ELSEIF(INDVAR==EP) THEN
                    REFLEN=RADU2*2.
                    EPINL=MAX(0.164*(ENINL**1.5)/(0.1*REFLEN),0.009)
                    CALL FN1(VAL,EPINL)
                  ENDIF
                ENDIF
              ELSEIF(LSYS) THEN
C ... Settings for System curve calculations
C ... Calculate Patch area at the beginning only
                IF(ISTEP==1.AND.ISWEEP==FSWEEP.AND.
     1            .NOT.LAREA(IC).AND.IZ>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
                L0VAL=L0F(VAL); L0DEN1=L0F(DEN1)
                L0P1=L0F(P1)
                IF(INDVAR==1) THEN
                  SETVAL=SYSVEL*IFUNC
C ... Set flow-direction velocity
                ELSEIF(INDVAR==3.OR.INDVAR==5.OR.INDVAR==7) THEN
                  SETVAL=SYSVEL*IVELD
                ENDIF
                L0XG=L0F(XG2D); L0YG=L0F(YG2D); L0ZG=L0F(ZGNZ); IPW=0
                DO IX=IXF,IXL
                  DO IY=IYF,IYL
                    ICELL=IY+NY*(IX-1)
                    IF(MASKPT(IPAT)) THEN ! skip cells not covered by object
                      IPW=IPW+1
                      IF(.NOT.L2DVRO(IPW)) THEN
                        F(L0VAL+ICELL)=0.0; CYCLE
                      ENDIF
                    ENDIF
                    IF(IFUNC<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=='X') THEN
                      INDV=3
                    ELSEIF(INJDIR=='Y') THEN
                      INDV=5
                    ELSE
                      INDV=7
                    ENDIF
                    IF(INDVAR==1) THEN
                      F(L0VAL+ICELL)=SETVAL*F(L0DEN1+ICELL)
                    ELSE
                      IF(ICFAN/=2.OR.INDV==INDVAR)
     1                          F(L0VAL+ICELL)=SETVAL
                    ENDIF
                    IF(ICFAN==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=='HTBX'.AND.NPATCH(1:2)=='TM') THEN
            ILEN=0
            IF(ISINZZ(NPATCH(3:3))) ILEN=1
            IF(ILEN==1.AND.ISINZZ(NPATCH(4:4))) ILEN=2
            IF(ILEN==2.AND.ISINZZ(NPATCH(5:5))) ILEN=3
            IF(ILEN==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==17.AND.INDVAR<=W2) THEN
C ... Set VALues for Fan-Matching
          IF(LMCH.AND.IC>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),IFUPDT,VELFAN(IC))
         if(dbg) write(14,'('' isweep'',i5,'' iz'',i3,'' c,v,vf,dp,a '',
     1  1p5e12.3)') isweep,iz,coset(ic),velarr(ic),velfan(ic)/parea(ic),
     1                                    fandp(ic)/parea(ic),parea(ic)
          ENDIF
        ENDIF
        RETURN
C***************************************************************
C
 
      ELSEIF(IGR==19) THEN
C ... ------- Section 2 ----- Start of sweep
        IF(ISC==2) THEN
          IF(LMCH.AND.ISTEP==1.AND.ISWEEP==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==FSWEEP) THEN
              IF(.NOT.DONE) THEN
                IF(MYID==0) CALL WRIT40
     +                 ('System characteristic option activated  ')
                CALL WRITBL
                SYSAR=0.0; SYSVEL=0.0; DPFAN=1.0E-10; DPSYS=1.E-10
                PHIGH=1.E-10; PLOW=0.0
                ITRMX=LSWEEP/(INPT+2); ITRCUT=ITRMX
                FLRTFX=AMAX1(FLRMIN,1.0)
                FLRTUP=(FLRMAX-FLRMIN)/(INPT+1)
                IF(MYID==0) THEN
                  NMFIL(53)='hotdat.csv'; CALL OPENFL(53)
                  WRITE(LUNIT(53),205)
                ENDIF
                DONE=LSYS
                DONEWR=.NOT.DONE
                ISCAL(1)=1
                IV=0
              ENDIF
            ELSE
              IF(INDVAR==P1) THEN ! update system flow rate
                IF(MYID==0) WRITE(LUNIT(53),206) ISWEEP,DPSYS,FLRTFX
                 IF(ISWEEP==ITRCUT.OR.(TOTRES(1)FSWEEP+1.AND.ISWEEP>ITRCUT0+5)) THEN
                  IV=IV+1
                  ITRCUT0=ISWEEP
                  VFLRT(IV)=FLRTFX
                  SPDRP(IV,1)=DPSYS
                  IF(IV0.0) THEN
              SYSVEL=FLRTFX/(3600.0*SYSAR)
            ELSE
              SYSVEL=0.0
              IF(LPAR) IFUNC=2
            ENDIF
          ELSE
C ... Fan Matching
            IF(LMCH) THEN
C.... Only do matching on visit for Pressure
              IF(INDVAR>P1) RETURN
C ... Check that for unsteady fan area has been set
              IF(ISTEP>1) THEN
                DO 1923 IT=1,IFANT
                  PAREA(IT)=AMAX1(PAREA(IT),TINY)
 1923           CONTINUE
              ENDIF
              IF(.NOT.DONE) THEN
                IF(MYID==0) THEN
                  CALL WRITST
                  CALL WRIT40('Fan Matching option activated')
                ENDIF
                DONE=LMCH
              ENDIF
              if(dbg.and.isweep==fsweep) then
                open(201,file='fan.csv',status='unknown',iostat=ios)
                close(201,status='delete',iostat=ios)
                open(201,file='fan.csv',status='new',iostat=ios)
                rewind 201
                write(201,'(
     1'' iswp,   velfan,        velnom,        velval,        velarr,'',
     1''        pdrp,          slope,         flarr'')')
              endif
              IF((ISWEEP>1.OR.QEQ(FIINIT(P1),READFI)).AND.
     1            (MOD(ISWEEP,IFUPDT)==0.OR.ISWEEP==FSWEEP)) THEN
                IF(ISTEP==FSTEP.AND.ISWEEP==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
C... Set CO to FIXFLU on 1st sweep of restart and set VAL to pressure drop
                      COSET(IT)=FIXFLU
                      VELARR(IT)=0.0
 1924               CONTINUE
                    ZEROED=.TRUE.
                  ENDIF
                ELSE
                  DO 192 IT=1,IFANT
C...  for parallel sum over processors
                    IF(LPAR) THEN
                      VOLF0=FANVEL(IT)*PAREA(IT); RDUM(1)=VOLF0
                      FAREA=PAREA(IT);            RDUM(2)=FAREA
                      FLFAN=FLARR(IT);            RDUM(3)=FLFAN
                      SDP=SYSDP(IT);              RDUM(4)=SDP
                      CALL RGSUMN(RDUM,4) ! sum over all processors
                      FAREA=RDUM(2)
                      VOLF0=RDUM(1); VOLF0=VOLF0/FAREA
                      FLFAN=RDUM(3)
                      SDP=RDUM(4)
                    ELSE
                      VOLF0=FANVEL(IT)
                      FAREA=PAREA(IT)
                      FLFAN=FLARR(IT)
                      SDP=SYSDP(IT)
                    ENDIF
C...  use possibly-summed values
                    PDRP=SDP/(FAREA+TINY)
                    IF(PDRP>DPMAX(IT)) THEN
C....   reduce flow rate by 10% when DP too big
                      FLFAN=FLFAN*.9
                    ELSEIF(PDRP==0.0) THEN
                      FLFAN=PRINTR(2,1,IT,DPMAX(IT)/2.0)
                    ELSE
                      FLFAN=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=FLFAN/(3600.0*FAREA*GEOMF(IT))
                    FLRPL1=FLFAN*1.01
                    VELPL1=FLRPL1/(3600.0*FAREA*GEOMF(IT))
                    PRPL1=PRINTR(1,2,IT,FLRPL1)
                    FLRM1=FLFAN/1.01
                    VELM1=FLRM1/(3600.0*FAREA*GEOMF(IT))
                    PRM1=PRINTR(1,2,IT,FLRM1)
                    SLOPE=((PRM1-PRPL1)+TINY)/((VELPL1-VELM1)+TINY)
C... Internal fan, or rectangular fan at domain edge
                    COSET(IT)=SLOPE*10000.0
                    VELVAL=VELNOM+PDRP/SLOPE/10000.0
C.... apply linear relaxation to fan velocity change
                    VELARR(IT)=(1.0-RELX)*VOLF0+RELX*VELVAL
                    FLARR(IT)=FLFAN*PAREA(IT)/FAREA
                    if(dbg) write(201,'(i4,'','',7(1pe14.5,'',''))')
     1    isweep,volf0,velnom,velval,velarr(it),pdrp,slope,flarr(it)
  192             CONTINUE
                ENDIF
              ENDIF
  205         FORMAT(3X,'ISWP,',8X,'SDP,',7X,'HFLRT')
  206         FORMAT(1X,I6,',',1PE11.3,',',1X,1PE11.3)
            ENDIF
          ENDIF
C   * ------------------- SECTION 4 ---- Start of iterations over slab.
        ELSEIF(ISC==4) THEN ! add fan source to velocity residual normalisation
          IF(LMCH.AND.(IZ==NZ.AND.(INDVAR==U1.OR.INDVAR==V1)).OR.
     1                                  (IZ==NZ-1.AND.INDVAR==W1)) THEN
            DO IPT=1,NUMPAT
              DO IC=1,IFANT
                IF(NAMPAT(IPT)==NAMPCH(IC).OR.NAMPAT(IPT)(4:)==
     1                                                 NAMPCH(IC)) THEN
                  IF(INDVAR==U1.AND.NX>1) THEN
                    CALL GETSO(IPT,U1,U1SO)
                    IF(QNE(U1SO,-999.)) F(L0TTFLM+IST(U1))=
     1                                      F(L0TTFLM+IST(U1))+ABS(U1SO)
                  ENDIF
                  IF(INDVAR==V1.AND.NY>1) THEN
                    CALL GETSO(IPT,V1,V1SO)
                    IF(QNE(V1SO,-999.)) F(L0TTFLM+IST(V1))=
     1                                      F(L0TTFLM+IST(V1))+ABS(V1SO)
                  ENDIF
                  IF(INDVAR==W1.AND.NZ>1) THEN
                    CALL GETSO(IPT,W1,W1SO)
                    IF(QNE(W1SO,-999.)) F(L0TTFLM+IST(W1))=
     1                                      F(L0TTFLM+IST(W1))+ABS(W1SO)
                  ENDIF
                ENDIF
              ENDDO
            ENDDO
          ENDIF
C ... ------- Section 6 ----- Finish of iz slab
        ELSEIF(ISC==6) THEN
          IF(ISWEEP==FSWEEP.AND.ISTEP==FSTEP) 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==2) THEN
            CALL HILOXY(L0F(P1),P1MX,P1MN,ILMX,ILMN,.false.)
            PHIGH=AMAX1(PHIGH,P1MX)
            IF(IZ==1) PLOW=PHIGH
            PLOW=AMIN1(PLOW,P1MN)
            IF(LPAR.AND.IZ==NZ) THEN
              CALL PAR_MAXR(PHIGH); CALL PAR_MINR(PLOW)
            ENDIF
          ENDIF
          IF((ISWEEP==LSWEEP.OR.ENUFSW).AND.ITEM1>0) THEN
            IF(IZ==1) THIGH=-GREAT
            CALL HILOXY(L0F(ITEM1),TM1MX,TM1MN,ILMX,ILMN,.FALSE.)
            IF(TM1MX>=THIGH) THEN
              IYLCMX=MOD(ILMX-1,NY)+1
              IXLCMX=(ILMX-IYLCMX)/NY+1
              IZLCMX=IZ
              THIGH=TM1MX
            ENDIF
            IF(IZ==1) TLOW=THIGH
            TLOW=AMIN1(TLOW,TM1MN)
            IF(LPAR.AND.IZ==NZ) THEN
              IZM=IZLCMX; IXM=IXLCMX; IYM=IYLCMX
              IJKM= IGL_IJK(IXM,IYM,IZM,0,0,0) ! get global IJK of max TEM1
              ITMAX=0.0; ITMAX(MYID+1)=IJKM
              CALL IGSUMN(ITMAX,NPROC)
              TMAX=0.0; TMAX(MYID+1)=THIGH ! load TMAX with local Thigh
              CALL RGSUMN(TMAX,NPROC)      ! sum over all procs
              IPROCM=MAXLOC(TMAX)          ! get processor number with largest TEM1
              IJKM=ITMAX(IPROCM(1))        ! get global IJK of overall Thigh
              NXNYM=NX0*NY0; NYM=NY0
              IZLCMX=1+(IJKM-1)/NXNYM; IJ=1+MOD(IJKM-1,NXNYM) ! recover IX,IY,IZ
              IXLCMX=1+(IJ-1)/NYM; IYLCMX=1+MOD(IJKM-1,NYM)
              THIGH=TMAX(IPROCM(1))        ! put global TMAX into Thigh
              TMAX=0.0; TMAX(MYID+1)=TLOW  ! load TMAX with local Tlow
              CALL RGSUMN(TMAX,NPROC)      ! sum over all procs
              IPROCM=MINLOC(TMAX)          ! get processor number with smallest TEM1
              TLOW=TMAX(IPROCM(1))         ! put global TMIN into Tlow
            ENDIF
          ENDIF
C   * ------------------- SECTION 7 ---- Finish of sweep.
        ELSEIF(ISC==7) THEN
          IF(LMCH.OR.LSYS) THEN
C... Pressure drop through the system
            IF(LMCH) THEN
              DO 1974 IT=1,IFANT
C... store pressure drop and density for fan before zeroing
                IF(NEZ(FDENA(IT))) THEN
                  SYSDP(IT)=FANDP(IT)
                  FDENS(IT)=FDENA(IT)/NSLB(IT)
                  FANVEL(IT)=VELFAN(IT)/PAREA(IT)
                ENDIF
                IF(ISWEEP/=LSWEEP) THEN
                  FANDP(IT)=0.0
                  FDENA(IT)=0.0
                  VELFAN(IT)=0.0
                ENDIF
 1974         CONTINUE
            ELSE
              DPINT=PHIGH-PLOW
              DPFAN=AMAX1(ABS(DPINT),DPFAN)
              PHIGH=PLOW
              DPSYS=DPFAN
              DPFAN=0.0
            ENDIF
            LOG1=.FALSE.
            IF(STEADY) THEN
              LOG1=((MOD(ISWEEP,IFREQ)==0.AND.CSG1(1:2)=='SW') .OR.
     1           MOD(ISWEEP,NPRINT)==0) .OR.ENUFSW.OR.ISWEEP==LSWEEP
            ELSE
              LOG1= (MOD(ISTEP,IFREQ)==0 .OR.MOD(ISTEP,NTPRIN)==0) .AND.
     1          (ENUFSW.OR.ISWEEP==LSWEEP)
            ENDIF
            IF(LMCH.AND.LOG1.AND..NOT.NULLPR) THEN
C ... Overall system flow-rate (source of R1)
              FLRTS=0.0
C.... sum flow rates over all processors
              IF(LPAR) THEN
                DO IR=1,GD_NUMPAT     ! loop over global patches
                  IPR=GD_INDPAT(Ir,1) ! Local Index of Patch
                  IF(IPR>0) THEN      ! on this processor
                    CALL GETSO(IPR,R1,SORCE)
                  ELSE                ! not on this processor
                    SORCE=0.0
                  ENDIF
                  IF(SORCE>0.0) FLRTS=FLRTS+SORCE ! sum inflow
                ENDDO
                CALL  GLSUM(FLRTS) ! sum over all processors
              ELSE
                DO IR=1,NUMREG
                  CALL GETSO(IR,R1,SORCE)
                  IF(SORCE>0.0) FLRTS=FLRTS+SORCE ! sum inflow
                ENDDO
              ENDIF
              IF(MYID==0) THEN
                CALL WRITST
                CALL WRITBL
                IF(STEADY)THEN
                  WRITE(LUPR1,
     1           '('' Operating point data at sweep = '',I6)') ISWEEP
                ELSE
                  WRITE(LUPR1,
     1           '('' Operating point data at sweep = '',I6)') ISWEEP
                  WRITE(LUPR1,'('' of time step = '',I6)') ISTEP
                  WRITE(LUPR1,'('' at time = '',1PE12.3,'' s'')') TIM
                ENDIF
                CALL WRITBL
              ENDIF
              IF(LEZ(BUOYD)) BUOYD=1.189
C....  No System flow rate if no external fans
              IFUNC1=2
              DO IT=1,IFANT
                IF(IFNARR(IT,2)/=2) IFUNC1=0
              ENDDO
              IF(MYID==0) THEN
                IF(FLRTS>0.0) THEN
                  WRITE(LUPR1,'(1X,1A,1P1E12.4,1A,1P1E12.4,1A)')
     1      'System mass flow rate: ',FLRTS,' kg/s',FLRTS*3600.0,' 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
              ENDIF
              DO 1972 IT=1,IFANT
C.... for parallel, sum flow rates etc over all processors
                IF(LPAR) THEN
                  FLTOT=FLARR(IT);          RDUM(1)=FLTOT
                  FAREA=PAREA(IT);          RDUM(2)=FAREA
                  FAND=FDENS(IT)*PAREA(IT); RDUM(3)=FAND
                  SDP=SYSDP(IT);            RDUM(4)=SDP
                  CALL RGSUMN(RDUM,4) ! sum ovr processors
                  FLTOT=RDUM(1);FAREA=RDUM(2);FAND=RDUM(3);SDP=RDUM(4)
                  FAND=FAND/FAREA; SDP=SDP/FAREA
                ELSE
                  FLTOT=FLARR(IT)
                  FAREA=PAREA(IT)
                  FAND=FDENS(IT)
                  SDP=SYSDP(IT)
                  SDP=SDP/FAREA
                ENDIF
                IF(MYID==0) THEN
                  WRITE(LUPR1,'(1X,4A)')
     1               'Fan Name: ',NAMARR(IT),' Fan Type: ',TYPARR(IT)
                  FANVL=FLTOT/(3600.0*FAREA*GEOMF(IT))
                  FANFLR=FLTOT*FAND
                  WRITE(LUPR1,'(1X,2A,1P1E12.4,1A,1P1E12.4,1A)')
     1           NAMARR(IT),' mass flow rate:   ',FANFLR/3600.,' kg/s ',
     1                                             FANFLR,' kg/h'
                  WRITE(LUPR1,'(1X,2A,1P1E12.4,1A,1P1E12.4,1A)')
     1            NAMARR(IT),' volume flow rate: ',FLTOT/3600.,' m^3/s',
     1                                             FLTOT,' m^3/h'
                  WRITE(LUPR1,'(1X,2A,1P1E12.4,1A)')
     1             NAMARR(IT),' velocity:         ',ABS(FANVL),' m/s'
                  WRITE(LUPR1,'(1X,2A,1P1E12.4,1A)')
     1             NAMARR(IT),' pressure drop:    ',SDP,' Pa'
                  CALL WRITBL
                ENDIF
 1972         CONTINUE
              IF(.NOT.STEADY) THEN
                CALL WRITBL
              ENDIF
              CALL WRITST
            ENDIF
          ENDIF
          IF(ISWEEP==LSWEEP.OR.ENUFSW) THEN
C... print sweep
            IF(ITEM1>0.AND.MYID==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.AND.MYID==0) THEN
              CALL CLOSFL(52)
              NMFIL(52)='sysdat.csv'
              CALL OPENFL(52)
              LUG=LUNIT(52)
              IF(LHOTDA==1) THEN
                WRITE(LUG,'(5X,1A,10X,1A)')'SVOL,','SDP'
                LHOTDA=2
                DO IV=1,INPT+2
                  IF(VFLRT(IV)>0.0.AND.SPDRP(IV,1)>0.0) THEN
                    WRITE(LUG,'(1X,1PE12.3,'','',1PE12.3)') VFLRT(IV),
     1                                                    SPDRP(IV,1)
                  ENDIF
                ENDDO
                CALL CLOSFL(52)
                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
      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 'satear'
      INCLUDE 'pltcfile'
      COMMON /WORDI1/ NWDS,NCHARS(20),NSEMI,H,NLINES
      COMMON /WORDL1/ ERROR,MORWDS
      COMMON /WORDC1/ WD(20),INLINE
      INTEGER         H
      LOGICAL         ERROR,MORWDS,OPTZZZ
      CHARACTER       WD*20, INLINE*120
      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,LINE1
      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/=0) THEN
C If error, attempt to read default fandata file
        NMFIL(IREF)=NMSAV
        CALL OPENZZ(IREF,IERR)
        IF(IERR/=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); REWIND(LUFD)
   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)==0) GOTO 10
      READ(LUFD,*) INUM
      PRDATA(ITYP,1,1)=INUM; PRDATA(ITYP,2,1)=INUM
      I=0
      DO WHILE(.TRUE.)
        READ(LUFD,'(A)',IOSTAT=IOS) LINE1
        IF(IOS/=0) EXIT
        ICOM=INDEX(LINE1,'!'); IF(ICOM>0) LINE1(ICOM:)=' ' ! allow ! as comment
        IF(LINE1==' '.OR.LINE1(1:1)=='*'.OR.LINE1(1:1)=='#') CYCLE ! skip blank lines & lines starti
        I=I+1; IF(I>INUM) EXIT
        LL=LENGZZ(LINE1)
        CALL SPLTZZZ(LINE1(1:LL),WD,NWDS,NCHARS,LL,' ,;'//CHAR(9),20) ! split into words allowing ta
        PRDATA(ITYP,1,I+1)=RRDZZZ(1); PRDATA(ITYP,2,I+1)=RRDZZZ(2)
      ENDDO
  100 IF(INDEX(LINE,NAMFAN)==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 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE '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, SLD, ISACTIVECELLS
C ...
      NAMFUN='FNAREA'
C ... Determine type for Area calculation
      IF(ITYP==2) THEN ! East
        IF(BFC) THEN
          L0AR=L0B(5)
        ELSE
          L0AR=L0F(AEAST)
        ENDIF
        L0VEL=L0F(U1); IPLUS=NY
        IF(IXFF==NX) L0VAL=L0F(WEST(U1))
      ELSEIF(ITYP==3) THEN ! West
        IF(BFC) THEN
          L0AR=L0B(6)
        ELSE
          L0AR=L0F(AEAST)
        ENDIF
        L0VEL=L0F(WEST(U1)); IPLUS=-NY
        IF(IXFF==1) L0VEL=L0F(U1)
      ELSEIF(ITYP==4) THEN ! North
        IF(BFC) THEN
          L0AR=L0B(7)
        ELSE
          L0AR=L0F(ANORTH)
        ENDIF
        L0VEL=L0F(V1); IPLUS=1
        IF(IYFF==NY) L0VAL=L0F(SOUTH(V1))
      ELSEIF(ITYP==5) THEN ! South
        IF(BFC) THEN
          L0AR=L0B(8)
        ELSE
          L0AR=L0F(ANORTH)
        ENDIF
        L0VEL=L0F(SOUTH(V1)); IPLUS=-1
        IF(IYFF==1) L0VEL=L0F(V1)
      ELSEIF(ITYP==6) THEN ! High
        IF(BFC) THEN
          L0AR=L0B(9)
        ELSE
          L0AR=L0F(AHIGH)
        ENDIF
        L0VEL=L0F(W1); IPLUS=NX*NY
        IF(IZ==NZ) L0VEL=L0F(LOW(W1))
      ELSEIF(ITYP==7) THEN ! Low
        IF(BFC) THEN
          L0AR=L0B(10)
        ELSE
          L0AR=L0F(AHIGH)
        ENDIF
        L0VEL=L0F(LOW(W1)); IPLUS=-NX*NY
        IF(IZ==1) L0VEL=L0F(W1)
      ENDIF
      AREA=0.0
      L0XG=L0F(XG2D); L0YG=L0F(YG2D); L0ZG=L0F(ZGNZ)
      DENA=0.0
      TVOL=0.0
      IPW=0
      L0DEN1=L0F(DEN1)
      DO IX=IXFF,IXLL
        DO IY=IYFF,IYLL
          IF(MASKPT(IPAT)) THEN ! skip cells not covered by object
            IPW=IPW+1
            IF(.NOT.L2DVRO(IPW)) CYCLE
          ENDIF
          IF(.NOT.ISACTIVECELLS(IX,IY,IZ)) CYCLE
          I=IY+NY*(IX-1)
          IF(SLD(I)) CYCLE
          IF((IPLUS== 1   .AND.IY1).OR.
     1       (IPLUS==NY   .AND.IX1).OR.
     1       (IPLUS==NX*NY.AND.IZ1)) THEN
            IF(SLD(I+IPLUS)) CYCLE
          ENDIF
            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
        ENDDO
      ENDDO
      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,IFUPDT,VEL)
C ... This subroutine calculates and sets value for each Fan type
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE '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, MASKPT, L2DVRO,
     1        ISACTIVECELLS
      SAVE LBNMUC, LBNMVC, LBNMWC
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
c
      NAMFUN='FANVAL'
      DONE1 = ISWEEP > FSWEEP
      LSETV = .NOT.BFC.OR.IFNC==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==FSWEEP.AND..NOT.LDONE).OR.ISWEEP1.1E-10) VEL=VEL+ABS(VELO)*AREA
          PDENA=PDENA+DENA
          LDONE=.FALSE.
          ICHK=IR
        ENDIF
      ENDIF
C ... Fan Setting part
      L0P1=L0F(P1); L0VAL=L0F(VAL)
      L0NGBR=L0F(P1)
      IF(ITYP==7) THEN
        L0P1=L0F(HIGH(P1)); L0AREA=L0F(AHIGH)
      ELSEIF(ITYP==3) THEN
        L0P1=L0F(EAST(P1)); L0AREA=L0F(AEAST)
      ELSEIF(ITYP==5) THEN
        L0P1=L0F(NORTH(P1)); L0AREA=L0F(ANORTH)
        ELSEIF(ITYP==2) THEN
          L0NGBR=L0F(EAST(P1)); L0AREA=L0F(AEAST)
        ELSEIF(ITYP==4) THEN
          L0NGBR=L0F(NORTH(P1)); L0AREA=L0F(ANORTH)
        ELSEIF(ITYP==6) THEN
          L0NGBR=L0F(HIGH(P1)); L0AREA=L0F(AHIGH)
        ENDIF
C ... Set flow-direction velocity
 
      IF(IVAR>=3.OR.IVAR<=7) SETVAL=ABS(FVEL)*IDIRV
C ... BFC part
      IF(.NOT.LSETV) THEN
        UCRT=0.0; VCRT=0.0; WCRT=0.0
        IUVW=3*((IVAR-U1)/2)
        IF(IVAR==3) THEN
          UCRT=SETVAL
        ELSEIF(IVAR==5) THEN
          VCRT=SETVAL
        ELSEIF(IVAR==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
      L0XG=L0F(XG2D); L0YG=L0F(YG2D); L0ZG=L0F(ZGNZ); IPW=0
      DO IX=IXFF,IXLL
        DO IY=IYFF,IYLL
          ICELL=IY+NY*(IX-1)
          IF(LSETV) F(L0VAL+ICELL)=0.0
          IF(MASKPT(IPAT)) THEN ! skip cells not covered by object
            IPW=IPW+1
            IF(.NOT.L2DVRO(IPW)) CYCLE
          ENDIF
          IF(.NOT.ISACTIVECELLS(IX,IY,IZ)) CYCLE
          IPCL=ICELL
C ... Fan pressure drop different for internal and external
          IF(IFNC==2) THEN
C ... Internal
            DP=F(L0P1+ICELL)-F(L0NGBR+ICELL)
          ELSE
C ... External
            IF(ITYP==3) IPCL=IY+NY*IX
            IF(ITYP==5) IPCL=IY+1+NY*(IX-1)
            DP=F(L0P1+IPCL)
          ENDIF
          FDP=FDP+ABS(DP)*F(L0AREA+ICELL)
          IF(ISWEEP==FSWEEP.AND.ISTEP==FSTEP.AND.FIINIT(P1)==READFI)THEN
             F(L0VAL+ICELL)=-DP/FIXFLU
          ELSE
            IF(LSETV) F(L0VAL+ICELL)=SETVAL
          ENDIF
        ENDDO
      ENDDO
      END
c