c
C.... file name GXROTASO.HTM 02.01.00
C**** SUBROUTINE GXROTA is called from group 13 of GREX3, by setting
C the value ascribed to GROUND in the COVAL statement. The
C subroutine is entered only when the patch name begins with the
C character 'ROTA'.
C For BFC=T, the subroutine BFCROT is called from this subroutine.
C
C.... The following library cases make use of it:
C B524,744-752(bfctst)
C
SUBROUTINE GXROTA(BFC,NX)
INCLUDE 'farray'
INCLUDE 'cmnrot'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
COMMON/RDATA/RDATA1(17),RINNER,RDATA2(67)
COMMON/IDATA/NXX,NY,NZ,IDT1(117)
COMMON/F0/KXDX,KXXG,KXXU,KXDXG,KXS,KYDY,KYYG,KYR,KYR2,KYRV,KYRV2,
1IF01(4),KYDYG,IF02(9),KZDZG,IF03(3),KXYDX,KXYDY,KTZPRF,KXYR,
1IF04(6),KXYDXG,KXYDYG,IF05(61),KD7,KD8,KD9,KD10,IF06(59),KXYDZ,
1IF08(138)
DIMENSION PA(3),PB(3)
COMMON /LROT/LROTOR
LOGICAL BFC,LROTOR,DOIT
COMMON /UVWCOL/ IUC1,IVC1,IFI(7)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
c
C.... functions called in this subroutine are:
c fn2(y,x,a,b) y = a+b*x
c fn10(y,x1,x2,a,b1,b2) y = a+b1*x1+b2*x2
c fn34(y,x,a) y = y+a*x
c
C.... IROTAA is a switch for activating a reduced-pressure
C system. IROTAA=1 activates this option. IROTAA=0 by default
C so that the centripetal pressure is included.
c
C.... ANGVEL is angular speed, clockwise about axis PA -> PB.
c
NAMSUB = 'GXROTA'
IF(ISC.GE.12) THEN
IF(BFC) THEN
C----------------------------------------------------------BFC=T
CALL VECTOR(PA,ROTAXA,ROTAYA,ROTAZA)
CALL VECTOR(PB,ROTAXB,ROTAYB,ROTAZB)
CALL BFCROT(VAL,PA,PB,ANGVEL,IROTAA)
ELSE
C----------------------------------------------------------BFC=F
IF(INDVAR.EQ.U1 .OR. INDVAR.EQ.U2) THEN
CALL FNAV(VAL,INDVAR+2,'SOUTH')
IF(NX.GT.1) CALL FNAV(VAL,VAL,'EAST')
IF(.NOT.LROTOR) THEN
CALL FN2(VAL,VAL,0.,-2.0*ANGVEL)
ELSE
L0VAL=L0F(VAL); IZADD=(IZ-1)*NX*NY
DO IX=IXF,IXL
DO IY=IYF,IYL
I=(IX-1)*NY+IY; IJK=I+IZADD; IROTC=INROT(IJK)
IF(IROTC>0) THEN
F(L0VAL+I)=-F(L0VAL+I)*2.*ANGV(IROTC)
ENDIF
ENDDO
ENDDO
ENDIF
ELSEIF(INDVAR.EQ.V1 .OR. INDVAR.EQ.V2) THEN
IF(.NOT.LROTOR) THEN
CALL FNAV(VAL,INDVAR-2,'NORTH')
ELSE
IF(LROTOR) THEN
L0VAL=L0F(VAL); L0VEL=L0F(INDVAR-2); IZADD=(IZ-1)*NX*NY
IYLL=MIN(IYL,NY-1)
DO IX=IXF,IXL
DO IY=IYF,IYLL
I=(IX-1)*NY+IY; IJK=I+IZADD
IROTC=INROT(IJK); IROTP=INROT(IJK+1)
IF(IROTC>0.AND.IROTP==0) THEN ! North face next to free space
F(L0VAL+I)=0.5*(F(L0VEL+I)+F(L0VEL+I+1)-
1 ANGV(IROTC)*F(KYR+IY+1))
ELSEIF(IROTC==0.AND.IROTP>0) THEN ! free space next to South face
F(L0VAL+I)=0.5*(F(L0VEL+I)+F(L0VEL+I+1)+
1 ANGV(IROTP)*F(KYR+IY+1))
ELSEIF(IROTC/=IROTP) THEN ! interface between adjacent rotors
F(L0VAL+I)=0.5*(F(L0VEL+I)+F(L0VEL+I+1)-
1 (ANGV(IROTC)-ANGV(IROTP))*F(KYR+IY+1))
ELSEIF(IROTC>0.AND.IROTP>0) THEN ! inside a rotor
F(L0VAL+I)=0.5*(F(L0VEL+I)+F(L0VEL+I+1))
ELSE ! in free space
F(L0VAL+I)=0.0
ENDIF
ENDDO
ENDDO
ELSE
CALL FNAV(VAL,INDVAR-2,'NORTH')
ENDIF
ENDIF
IF(NX.GT.1) CALL FNAV(VAL,VAL,'WEST')
IF(.NOT.LROTOR) THEN
CALL FN2(VAL,VAL,0.0,2.0*ANGVEL)
IF(IROTAA.EQ.0) CALL FN34(VAL,RV2D,ANGVEL*ANGVEL)
ELSE
IZADD=(IZ-1)*NX*NY; L0VAL=L0F(VAL); L0R=L0F(RV2D)
DO IX=IXF,IXL
DO IY=IYF,IYL
I=(IX-1)*NY+IY; IJK=I+IZADD; IROTC=INROT(IJK)
IF(IROTC>0) THEN
F(L0VAL+I)=F(L0VAL+I)*2.*ANGV(IROTC)
IF(IROTAA==0)THEN
F(L0VAL+I)=F(L0VAL+I)+ANGV(IROTC)*ANGV(IROTC)*
1 F(L0R+I)
ENDIF
ELSE
F(L0VAL+I)=0.0
ENDIF
ENDDO
ENDDO
ENDIF
ELSEIF(INDVAR.EQ.IUC1) THEN
C.... X-direction colocated component:
CALL FN0(VAL,IVC1)
CALL FN2(VAL,VAL,0.,-2.0*ANGVEL)
ELSEIF(INDVAR.EQ.IVC1) THEN
C.... Y-direction colocated component:
CALL FN0(VAL,IUC1)
CALL FN2(VAL,VAL,0.0,2.0*ANGVEL)
IF(IROTAA.EQ.0) CALL FN34(VAL,RV2D,ANGVEL*ANGVEL)
ENDIF
ENDIF
ENDIF
NAMSUB = 'gxrota'
END