c

c.... File Name .....   gxswfan.htm ...  221201
C ... Swirl functions for Fans
      SUBROUTINE SWIRFAN
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
      PARAMETER (NPRP=50)
      COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15)
     1           ,ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY
     1           ,IRADZ,IVFOL
      COMMON/DRHODP/ITEMP,IDEN/DVMOD/IDVCGR
      COMMON/NAMFN/NAMFUN,NAMSUB
      COMMON/SWRFAN/RADU1,RADU2,GX0,GY0,GZ0,SWN,ICFAN,IFANN(NPRP)
     1             ,INJDIR
      COMMON/FANMS1/IFANT,VELARR(NPRP),NAMPCH(NPRP)
      CHARACTER*6 NAMFUN,NAMSUB,NAMPCH*16
      CHARACTER*14 COBJTP,OBNAME,CTEMP*8,INJDIR*1
      LOGICAL MASKPT
C
C***********************************************************************
      NAMSUB='SWRFAN'
      IXL=IABS(IXL)
C*****************************************************************
C--- GROUP 13. Boundary conditions and special sources
C                                       Index for Coefficient - CO
C                                       Index for Value       - VAL
      IF(IGR.EQ.13) THEN
        IF(ISC.EQ.12.OR.ISC.EQ.17) THEN
C------------------- SECTION 12 ------------------- value = GRND
C
C  USWIRL       ----- swirl velocity  & direction (+ve for clockwise
C                     into the solution domain)
C  GXO,GY0,GZO  ----- coordinates of burner axis
C  INJDIR       ----- direction of injection
C
          IF(MASKPT(IPAT)) THEN
            OBNAME=COBJTP(IPAT)
            IF(OBNAME(1:3).EQ.'FAN') THEN
C...  Decomposition of Swirl Velocity Vector
              CTEMP='^'//NPATCH
              CALL GETSDI(CTEMP,'ICFAN',ICFAN)
              CALL GETSDR(CTEMP,'XVEL1',XVEL1)
              CALL GETSDR(CTEMP,'YVEL1',YVEL1)
              CALL GETSDR(CTEMP,'ZVEL1',ZVEL1)
              CALL GETSDR(CTEMP,'FVEL1',FVEL1)
              CALL GETSDR(CTEMP,'DIAM1',DIAM1)
              CALL GETSDR(CTEMP,'DIAM2',DIAM2)
              CALL GETSDR(CTEMP,'GX0',GX0)
              CALL GETSDR(CTEMP,'GY0',GY0)
              CALL GETSDR(CTEMP,'GZ0',GZ0)
              CALL GETSDR(CTEMP,'SWIRLN',SWN)
              CALL GETSDC(CTEMP,'INJDIR',INJDIR)
              RADU1=DIAM1/2.0
              RADU2=DIAM2/2.0
              USWIRL=ABS(FVEL1)*SWN
              DO IC=1,IFANT
                IF(NPATCH.EQ.NAMPCH(IC)) THEN
                  USWIRL=VELARR(IC)*SWN
                  IFANN(IC)=ICFAN
                ENDIF
              ENDDO
            ENDIF
          ENDIF
        ENDIF
C ... If not swirl fan, return
        IF(ICFAN.NE.2) RETURN
        IF(ISC.EQ.12) THEN
C------------------- SECTION 12 ------------------- value = GRND
C
          IF(MASKPT(IPAT)) THEN
            OBNAME=COBJTP(IPAT)
            IF(OBNAME(1:3).EQ.'FAN') THEN
C
C...  X-direction
C
              IF(INJDIR.EQ.'X') THEN
C...  Axial injection in x direction through y-z plane
                IF(INDVAR.EQ.V1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0ZGG,L0F(ZGNZ)
     1                     ,L0YG,L0F(YG2D))
                  L0ZG=L0VAL
                  DO 13801 IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO 13801 IY=IYF,IYL
                      I=IY+IADD
                      F(L0ZG+I)=F(L0ZGG+IZ)
13801             CONTINUE
                  CALL GXSINA(L0VAL,L0YG,L0ZG,GY0,GZ0,NY)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0ZGG+IZ)-GZ0)**2+(F(L0YG+I)-GY0)**2
                      DIS=SQRT(DIS)
c ...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=F(L0VAL+I)*(-USWIRL)
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ELSEIF(INDVAR.EQ.W1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0ZGG,L0F(ZGNZ)
     1                     ,L0YG,L0F(YG2D))
                  L0ZG=L0VAL
                  DO 13802 IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO 13802 IY=IYF,IYL
                      I=IY+IADD
                      F(L0ZG+I)=F(L0ZGG+IZ)
13802             CONTINUE
                  CALL GXCOSA(L0VAL,L0YG,L0ZG,GY0,GZ0,NY)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0ZGG+IZ)-GZ0)**2+(F(L0YG+I)-GY0)**2
                      DIS=SQRT(DIS)
c...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=F(L0VAL+I)*USWIRL
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
c... P1 fixflu
                ELSEIF(INDVAR.EQ.P1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0ZGG,L0F(ZGNZ)
     1                     ,L0YG,L0F(YG2D))
                  L0DEN1=L0F(DEN1)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0ZGG+IZ)-GZ0)**2+(F(L0YG+I)-GY0)**2
                      DIS=SQRT(DIS)
c...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=FVEL1*F(L0DEN1+I)
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ELSEIF(INDVAR.EQ.U1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0ZGG,L0F(ZGNZ)
     1                     ,L0YG,L0F(YG2D))
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0ZGG+IZ)-GZ0)**2+(F(L0YG+I)-GY0)**2
                      DIS=SQRT(DIS)
c...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=XVEL1
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ENDIF
C
C ... Y-direction
C
              ELSEIF(INJDIR.EQ.'Y') THEN
C ... Axial injection in y direction through x-z plane
                IF(INDVAR.EQ.U1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0ZGG,L0F(ZGNZ)
     1                     ,L0XG,L0F(XG2D))
                  L0ZG=L0VAL
                  DO 13803 IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO 13803 IY=IYF,IYL
                      I=IY+IADD
                      F(L0ZG+I)=F(L0ZGG+IZ)
13803             CONTINUE
                  CALL GXSINA(L0VAL,L0XG,L0ZG,GX0,GZ0,NY)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0ZGG+IZ)-GZ0)**2+(F(L0XG+I)-GX0)**2
                      DIS=SQRT(DIS)
c ... Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=F(L0VAL+I)*USWIRL
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ELSEIF(INDVAR.EQ.W1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0ZGG,L0F(ZGNZ)
     1                     ,L0XG,L0F(XG2D))
                  L0ZG=L0VAL
                  DO 13804 IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO 13804 IY=IYF,IYL
                      I=IY+IADD
                      F(L0ZG+I)=F(L0ZGG+IZ)
13804             CONTINUE
                  CALL GXCOSA(L0VAL,L0XG,L0ZG,GX0,GZ0,NY)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0ZGG+IZ)-GZ0)**2+(F(L0XG+I)-GX0)**2
                      DIS=SQRT(DIS)
c ...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=F(L0VAL+I)*(-USWIRL)
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
c ... P1 fixflu
                ELSEIF(INDVAR.EQ.P1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0ZGG,L0F(ZGNZ)
     1                     ,L0XG,L0F(XG2D))
                  L0DEN1=L0F(DEN1)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0ZGG+IZ)-GZ0)**2+(F(L0XG+I)-GX0)**2
                      DIS=SQRT(DIS)
c ...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=FVEL1*F(L0DEN1+I)
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ELSEIF(INDVAR.EQ.V1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0ZGG,L0F(ZGNZ)
     1                     ,L0XG,L0F(XG2D))
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0ZGG+IZ)-GZ0)**2+(F(L0XG+I)-GX0)**2
                      DIS=SQRT(DIS)
c ...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=YVEL1
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ENDIF
c
c...  Z-direction
c
              ELSEIF(INJDIR.EQ.'Z') THEN
C...  Axial injection in z direction through y-x plane
                IF(INDVAR.EQ.U1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0YG,L0F(YG2D)
     1                     ,L0XG,L0F(XG2D))
                  CALL GXSINA(L0VAL,L0XG,L0YG,GX0,GY0,NY)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0XG+I)-GX0)**2+(F(L0YG+I)-GY0)**2
                      DIS=SQRT(DIS)
c ...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=F(L0VAL+I)*(-USWIRL)
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ELSEIF(INDVAR.EQ.V1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0YG,L0F(YG2D)
     1                     ,L0XG,L0F(XG2D))
                  CALL GXCOSA(L0VAL,L0XG,L0YG,GX0,GY0,NY)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0XG+I)-GX0)**2+(F(L0YG+I)-GY0)**2
                      DIS=SQRT(DIS)
c ...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=F(L0VAL+I)*USWIRL
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
c ... P1 fixflu
                ELSEIF(INDVAR.EQ.P1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0YG,L0F(YG2D)
     1                     ,L0XG,L0F(XG2D))
                  L0DEN1=L0F(DEN1)
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0XG+I)-GX0)**2+(F(L0YG+I)-GY0)**2
                      DIS=SQRT(DIS)
c ...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=FVEL1*F(L0DEN1+I)
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ELSEIF(INDVAR.EQ.W1) THEN
                  CALL SUB3(L0VAL,L0F(VAL),L0YG,L0F(YG2D)
     1                     ,L0XG,L0F(XG2D))
                  DO IX=IXF,IXL
                    IADD=NY*(IX-1)
                    DO IY=IYF,IYL
                      I=IY+IADD
                      DIS=(F(L0XG+I)-GX0)**2+(F(L0YG+I)-GY0)**2
                      DIS=SQRT(DIS)
c ...  Block centre area
                      IF(DIS.GE.RADU1.AND.DIS.LE.RADU2) THEN
                        F(L0VAL+I)=ZVEL1
                      ELSE
                        F(L0VAL+I)=0.0
                      ENDIF
                    ENDDO
                  ENDDO
                ENDIF
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      NAMSUB='swrfan'
      END
c