c

C name htbxgr.for................................................081117
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'
      INCLUDE '/phoenics/d_includ/parear'
      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),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)
      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,L2DVRO
      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: 081117')
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'
C#### JCL 07.06.16 add turbulence intensity for edge matched fan
              TURBN(IQ1)=0.0
              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))
C#### JCL 05.01.07 upper-case fan type
              LL=LENGZZ(TYPARR(IQ1))
              CALL UCASZZ(TYPARR(IQ1),LL)
              CALL GETSDC(NAMFAN,'PATFAN',NAMPCH(IQ1))
              IF(NAMARR(IQ1).NE.' ') IFANT=IFANT+1
C#### JCL 07.06.16 add turbulence intensity for edge matched fan
              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)
          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)
C#### JCL 09.01.13 count slabs separately for facetted and non-factetted                  
C####                  NSLB(IC)=K2-K1+1
                  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.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
C#### JCL 09.01.13 count slabs in facetted object otherwise fan mass flow rate printed wrong
                        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(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
                    NSLB(IC)=K2-K1+1
                  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)
C#### JCL 07.06.16 next section only for swirl fan
              IF(LMCH.AND.ICFAN==2) THEN
C... Fixflu for P1 for circular Fans
                IF(INDVAR.EQ.P1) THEN
                  CALL SUB3(L0XG,L0F(XG2D), L0YG,L0F(YG2D)
     1                      ,L0ZG,L0F(ZGNZ))
                  L0DEN1=L0F(DEN1); L0VAL=L0F(VAL); IPW=0
                  DO IX=IXF,IXL
                    DO IY=IYF,IYL
                      I=IY+NY*(IX-1)
C#### JCL 07.06.16 use object mask instead of radius to match SWIRFAN
                      IF(MASKPT(IPAT)) THEN
                        IPW=IPW+1
                        IF(.NOT.L2DVRO(IPW)) THEN
                          F(L0VAL+I)=0.0; CYCLE
                        ENDIF
                      ENDIF  
C#### JCL 07.06.16 use object mask instead of radius to match SWIRFAN
C####c... dis, distance from circle centre to cell centre
C####                      IF(INJDIR.EQ.'X') THEN
C####                        DIS=(F(L0YG+I)-GY0)**2+(F(L0ZG+IZ)-GZ0)**2
C####                      ELSEIF(INJDIR.EQ.'Y') THEN
C####                        DIS=(F(L0XG+I)-GX0)**2+(F(L0ZG+IZ)-GZ0)**2
C####                      ELSE
C####                        DIS=(F(L0XG+I)-GX0)**2+(F(L0YG+I)-GY0)**2
C####                      ENDIF
C####                      DIS=SQRT(DIS)
C####                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=F(L0DEN1+I)*VELARR(IC)*IFUNC
C####                      ELSE
C####                        F(L0VAL+I)=0.0
C####                      ENDIF
                    ENDDO
                  ENDDO
C#### JCL 07.06.16 set inlet KE, EP values for matched fan at domain edge
                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.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
                L0XG=L0F(XG2D); L0YG=L0F(YG2D); L0ZG=L0F(ZGNZ); IPW=0
                DO IX=IXF,IXL
                  DO IY=IYF,IYL
                    ICELL=IY+NY*(IX-1)
C#### JCL 07.06.16 use object mask instead of radius to match SWIRFAN
                    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.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
C####                   !  DIS=(F(L0YG+ICELL)-GY0)**2+(F(L0ZG+IZ)-GZ0)**2
                    ELSEIF(INJDIR.EQ.'Y') THEN
                      INDV=5
C####                      DIS=(F(L0XG+ICELL)-GX0)**2+(F(L0ZG+IZ)-GZ0)**2
                    ELSE
                      INDV=7
C####                      DIS=(F(L0XG+ICELL)-GX0)**2+(F(L0YG+ICELL)-GY0)**2
                    ENDIF
C####                    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
C####                    IF((ICFAN.EQ.2).AND.
C####     1                 (DIS.LT.RADU1.OR.DIS.GT.RADU2))
C####     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
C#### JCL 05.01.07 add test for velocity to avoid conflict with
C####              T3 and THINPLT
        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; 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)
C#### JCL 12.12.16 change HOTDAT to HOTDAT.CSV
                NMFIL(53)='hotdat.csv'; 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
C#### JCL 12.12.16 TOTRES(1) condition should refer to RESFAC not 1.0, and 
C####              may well be met on first sweeps. 
C####              IF(ISWEEP.EQ.ITRCUT.OR.TOTRES(1).LT.1) THEN
!              IF(ISWEEP.EQ.ITRCUT) THEN
               IF(ISWEEP.EQ.ITRCUT.OR.(TOTRES(1).LT.RESFAC.AND.
     1                                            ISWEEP>FSWEEP+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
C#### JCL 05.01.07 set CO large to prevent movement on 1st restart sweep
                      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)
C#### JCL 07.06.16 use linearisation for all cases
C####                    IF(IFANN(it).NE.2.or.IFNARR(IT,2).EQ.2) THEN
C####C... Internal fan, or rectangular fan at domain edge
                    COSET(IT)=SLOPE*20.0
                    VELVAL=VELNOM+PDRP/SLOPE/20.0
                    VELARR(IT)=VELVAL
C####                    ELSE
C####                      COSET(IT)=ONLYMS
C####                      VELARR(IT)=VELARR(IT)-(VELNOM-VELARR(IT))*
C####     1                                                    (PDRP/SLOPE)
C####                    ENDIF
  192             CONTINUE
                ENDIF
              ENDIF
C#### JCL 12.12.16 change formats of HOTDAT to match CSV
C####  205         FORMAT(1X,4HISWP,9X,3HSDP,9X,4HFLRT)
C####  206         FORMAT(1X,I4,1X,1PE11.3,2X,1PE11.3)
  205         FORMAT(3X,'ISWP,',8X,'SDP,',7X,'HFLRT')
  206         FORMAT(1X,I6,',',1PE11.3,',',1X,1PE11.3)
            ENDIF
          ENDIF
C ... ------- Section 6 ----- Finish of iz slab
        ELSEIF(ISC.EQ.6) THEN
C#### JCL 12.02.10 check on FSWEEP+1 to allow restart with FSTEP>1
C####          IF(ISWEEP.LE.2) 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
C#### JCL 27.05.09 store pressure drop and density for fan before zeroing
                IF(NEZ(FDENA(IT))) THEN
                  SYSDP(IT)=FANDP(IT)
                  FDENS(IT)=FDENA(IT)
                ENDIF
C#### JCL 10.09.13 do need to zero when ENUFSW=T on early termination
C####                IF(ISWEEP.NE.LSWEEP.AND..NOT.ENUFSW) THEN
                IF(ISWEEP.NE.LSWEEP) 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))
C#### JCL 27.05.09 use stored value of density
C####                FANFLR=FLARR(IT)*FDENA(IT)/NSLB(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)')
C#### JCL 27.05.09 use stored value of pressure drop
C####     1             NAMARR(IT),' pressure drop:    ',FANDP(IT)
     1             NAMARR(IT),' pressure drop:    ',SYSDP(IT)
     1            ,' Pa'
                CALL WRITBL
 1972         CONTINUE
              IF(.NOT.STEADY) THEN
                CALL WRITBL
              ENDIF
              CALL WRITST
            ENDIF
          ENDIF
C####          IF(ISWEEP.EQ.LSWEEP.OR.ENUFSW) THEN
          IF(ISWEEP.EQ.LSWEEP) THEN
C... print sweep
            IF(ITEM1.GT.0.AND.MYID.EQ.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)
C#### JCL 07.06.16 change name to lower-case and add .csv. 
              NMFIL(52)='sysdat.csv'
              CALL OPENFL(52)
              LUG=LUNIT(52)
C#### JCL 12.12.16 add mssing ',' after SVOL
              IF(LHOTDA.EQ.1)WRITE(LUG,'(5X,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,1PE12.3,'','',1PE12.3)') VFLRT(IV),
     1                                                    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')
C#### gin 30.03.05 New error handling
      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'
      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.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)
C#### JCL 08.11.17 LUDF on next line should be LUFD
      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).EQ.0) GOTO 10
      READ(LUFD,*) INUM
      PRDATA(ITYP,1,1)=INUM; PRDATA(ITYP,2,1)=INUM
C#### JCL 23.02.17 allow for space,; and tab as delimiters in FANDATA file      
C####      DO 20 I=1,INUM
C####         READ(LUFD,*) RVL,RDP
C####         PRDATA(ITYP,1,I+1)=RVL
C####         PRDATA(ITYP,2,I+1)=RDP
C####   20 CONTINUE
      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 starting #
        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 tab 
        PRDATA(ITYP,1,I+1)=RRDZZZ(1); PRDATA(ITYP,2,I+1)=RRDZZZ(2)
      ENDDO 
  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***********************************************************************
C#### JCL 05.01.07 return average velocity
      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, SLD
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 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
          I=IY+NY*(IX-1)
C#### JCL 07.06.16 only use object mask, and proper test on solid
C####          IF(ICFAN.EQ.2) THEN
C####            IF(INJDIR.EQ.'X') THEN
C####              DIS=(F(L0YG+I)-GY0)**2+(F(L0ZG+IZ)-GZ0)**2
C####            ELSEIF(INJDIR.EQ.'Y') THEN
C####              DIS=(F(L0XG+I)-GX0)**2+(F(L0ZG+IZ)-GZ0)**2
C####            ELSE
C####              DIS=(F(L0XG+I)-GX0)**2+(F(L0YG+I)-GY0)**2
C####            ENDIF
C####            DIS=SQRT(DIS)
C####            IF(DIS.LT.RADU1.OR.DIS.GT.RADU2) GOTO 1332
C####          ENDIF
C####C#### JCL 08.01.07 put limit on dens to exclude solid cells
C####          IF(F(L0DEN1+I).LE.2.0) THEN
          IF(.NOT.SLD(I)) 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
        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)
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, MASKPT, L2DVRO
      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
C#### JCL 05.01.07 return average velocity in VEL
          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#### JCL 05.01.07 do this each sweep in case run finishes early
C####C ... Calculate Product of Patch area and Density at the end only
C####      IF(ISWEEP.EQ.LSWEEP.AND.LDONE) THEN
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

C####      IF(EQZ(FVEL)) FVEL=0.01
      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
      L0XG=L0F(XG2D); L0YG=L0F(YG2D); L0ZG=L0F(ZGNZ); IPW=0
      DO IX=IXFF,IXLL
        DO IY=IYFF,IYLL
          ICELL=IY+NY*(IX-1)
C#### JCL 07.06.16 use object mask
          IF(MASKPT(IPAT)) THEN ! skip cells not covered by object
            IPW=IPW+1
            IF(.NOT.L2DVRO(IPW)) THEN
              IF(LSETV) F(L0VAL+ICELL)=0.0; CYCLE
            ENDIF
          ENDIF  
          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
C#### JCL 07.06.16 use object mask
C####          IF(ICFAN.EQ.2) THEN
C####          IF(INJDIR.EQ.'X') THEN
C####            DIS=(F(L0YG+ICELL)-GY0)**2+(F(L0ZG+IZ)-GZ0)**2
C####          ELSEIF(INJDIR.EQ.'Y') THEN
C####            DIS=(F(L0XG+ICELL)-GX0)**2+(F(L0ZG+IZ)-GZ0)**2
C####          ELSE
C####            DIS=(F(L0XG+ICELL)-GX0)**2+(F(L0YG+ICELL)-GY0)**2
C####          ENDIF
C####          DIS=SQRT(DIS)
C####          IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
C####            IF(LSETV) F(L0VAL+ICELL)=SETVAL
C####          ELSE
C####            IF(LSETV) F(L0VAL+ICELL)=0.0
C####          ENDIF
C####          ELSE
          IF(LSETV) F(L0VAL+ICELL)=SETVAL
C####          ENDIF
        ENDDO
      ENDDO  
      END
c