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