c

C.... file name GXMISCSO.F                                 261109
C**** SUBROUTINE GXPOLR is called from group 13 of GREX3, by setting
C     the value ascribed to GROUND in the COVAL statement, and is
C     entered only when the patch name with the 2nd to 4th
C     characters of 'POL'.
C
C.... The dummy NP is the first two characters of the patch name
C     used in the SATELLITE; UP should be used for x-direction momentum
C     equations, and VP for the y-direction momentum equations.
C
C.... The library cases 222-226,237 make use of it.
C
      SUBROUTINE GXPOLR(NP)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON/IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
     1       ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
      INTEGER VAL,CO,WALDIS,PATGEO
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
C
c     fn25(y,a)                        y = a*y
c     fn41(y,x,a,b)                    y = sin(a*x+b)
c     fn42(y,x,a,b)                    y = cos(a*x+b)
c     POLRA is the flow speed.
c
      NAMSUB = 'GXPOLR'
c------------------------------ val = grnd1 or setspeed
      IF(ISC.EQ.13) THEN
        IF(NP.EQ.'UP'.AND.(INDVAR.EQ.U1.OR.INDVAR.EQ.U2)) THEN
          CALL FN41(VAL,XU2D,1.0,0.0)
          CALL FN25(VAL,-POLRA)
        ELSE
          IF(NP.EQ.'VP'.AND.(INDVAR.EQ.V1.OR.INDVAR.EQ.V2)) THEN
            CALL FN42(VAL,XG2D,1.0,0.0)
            CALL FN25(VAL,POLRA)
          ENDIF
        ENDIF
      ENDIF
      NAMSUB = 'gxpolr'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GXPDRP is called from group 13 of EARTH, by setting the coefficient
C     to GRND-value in the SPEDAT command set for a porous 'plate' or
C     'thin-plate' object:
C                  SPEDAT(SET,PDROP_LAW,NamPat,R,GRND*),
C     where NamPat is the name of the object. This ascribes a pressure
C     drop formula, determined by the value of GRND*. as follows:
C       GRND1 - Dp = Const * U**Power
C       GRND2 - Dp = Const * 0.5 * Dens * U**2,
C     here U is the device or superficial velocity component across the plate.
C     The choice between superficial and device is made by a SPEDAT command:
C                  SPEDAT(SET,PDROP_VEL,NamPat,R,value)
C     where value = 0 for device velocity or = 1 for superficial velocity.
C     The coefficient for GRND1 and GRND2 is set by another SPEDAT:
C                  SPEDAT(SET,PDROP_COE,NamPat,R,value)
C     The power for GRND1 is set by the SPEDAT:
C                  SPEDAT(SET,PDROP_POW,NamPat,R,value)
C     Note, that this command will have an effect only if the porosity
C     of a 'thin-plate' object is set to a non-zero value by the command
C                  SPEDAT(SET,POROSITY,NamPat,R,Value),
C     where 0.< Value <1. specifies the plate active area.
C
C     A 'thin-plate' allows heat transfer through the plate. The nominal
C     thickness and material number are set by:
C                  COVAL(NamPat,PRPS, Material, Thickness)
C
C     A 'plate' does not allow heat transfer, but does allow different
C     heat sources on either side.
C
C     NOTE: If this routine is to be used in a PIL-based Q1, the following
C     rules should be observed:
C
C     A THIN-PLATE is defined by a PATCH with name starting PLT*.
C     NamPat in the SPEDATs above must be the same as the patch name.
C
C     A PLATE is defined by a PATCH with name starting PPD*, e.g. PPD*abcd
C     NamPat in the SPEDATs above must start with PLT*, e.g. PLT*abcd.
 
      SUBROUTINE GXPDRP
      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'
      COMMON/NAMFN/NAMFUN,NAMSUB
      LOGICAL FLUID1,MASKPT
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB= 'GXPDRP'
      IF(ISC.EQ.2 .OR. ISC.EQ.3) THEN
        L0CO = L0F(CO)
        L0VEL= L0F(INDVAR)
        CONST= GETPCN(0,IPAT,MASKPT(IPAT))
        IF(ISC.EQ.2) THEN
C.... GRND1 coefficient gives pressure drop proportional to the power of
C     the superficial velocity across the plate:
          POWER= GETPCN(1,IPAT,MASKPT(IPAT)) - 1.
          DO IX= IXF,IXL
            IAD= (IX-1)*NY
            DO IY= IYF,IYL
              I= IY + IAD
              VELPLT= F(L0VEL+I)+TINY
              F(L0CO+I)= CONST*ABS(VELPLT)**POWER
            ENDDO
          ENDDO
        ELSE
C.... GRND2 coefficient gives pressure drop proportional to the dynamic
C     pressure across the plate:
          LRHO12= ITWO(LRHO1,LRHO2,FLUID1(INDVAR))
          L0RHO = L0F(LRHO12)
          IF(INDVAR.EQ.U1.OR.INDVAR.EQ.U2) THEN
            LRHON=EAST(LRHO12)
          ELSEIF(INDVAR.EQ.V1.OR.INDVAR.EQ.V2) THEN
            LRHON=NORTH(LRHO12)
          ELSEIF(INDVAR.EQ.W1.OR.INDVAR.EQ.W2) THEN
            LRHON=ITWO(LRHO1H,LRHO2H,FLUID1(INDVAR))
          ENDIF
          L0RHON=L0F(LRHON)
          DO IX= IXF,IXL
            IAD= (IX-1)*NY
            DO IY= IYF,IYL
              I= IY + IAD
              IF(F(L0VEL+I).GT.0.0) THEN
                DENS=F(L0RHO+I)
              ELSE
                DENS=F(L0RHON+I)
              ENDIF
              F(L0CO+I)= 0.5*CONST*DENS*(ABS(F(L0VEL+I))+TINY)
            ENDDO
          ENDDO
        ENDIF
      ENDIF
      NAMSUB= 'gxpdrp'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GXDENDIF is called from group 13 of EARTH, for patch names beginning
C     DENDIF. It sets VAL to for velocities across faces having densities
C     on each side which differ by than 1.0
C
      SUBROUTINE GXDENDIF
      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'
      COMMON/NAMFN/NAMFUN,NAMSUB
      LOGICAL FLUID1
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB= 'DENDIF'
      L0CO  = L0F(CO)
      L0VAL = L0F(VAL)
      LRHO12= ITWO(LRHO1,LRHO2,FLUID1(INDVAR))
      L0RHO = L0F(LRHO12)
      IF(INDVAR.EQ.U1.OR.INDVAR.EQ.U2) THEN
        LRHON=EAST(LRHO12)
      ELSEIF(INDVAR.EQ.V1.OR.INDVAR.EQ.V2) THEN
        LRHON=NORTH(LRHO12)
      ELSEIF(INDVAR.EQ.W1.OR.INDVAR.EQ.W2) THEN
        LRHON=ITWO(LRHO1H,LRHO2H,FLUID1(INDVAR))
      ENDIF
      L0RHON=L0F(LRHON)
      DO 20 IX= IXF,IXL
        IAD= (IX-1)*NY
        DO 20 IY= IYF,IYL
          I= IY + IAD
          DIFF=ABS(F(L0RHO+I)  - F(L0RHON+I))
          IF(DIFF.GT.1.0) THEN
            F(L0VAL+I)= 0.0
          ELSE
            F(L0CO+I) = 0.0
          ENDIF
   20 CONTINUE
      NAMSUB= 'dendif'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GETPCN retrieves value of either coefficient (NGO=0), or power
C     (NGO=1) specified for a porous 'thin-plate' object. GETPCN is
C     called from GXPDRP.
      FUNCTION GETPCN(NGO,IPAT,VRPAT)
      INCLUDE '/phoenics/d_includ/patcmn'
      LOGICAL VRPAT,LTMP,FIRST
      CHARACTER NPATCH*8, MODEL*10, CTMP
C.... Preliminaries. Conver the name of a pressure drop patch into
C     the name of the correponding 'thin-plate' patch
      IF(VRPAT) THEN
        NPATCH= '^PLT*'//NAMPAT(IPAT)(5:7)
      ELSE
        NPATCH= 'PLT*'//NAMPAT(IPAT)(5:8)
      ENDIF
      GETPCN= 0.0
      FIRST = .TRUE.
      IF(NGO.EQ.0) THEN
        MODEL= 'PDROP_COE'
      ELSE
        MODEL= 'PDROP_PWR'
      ENDIF
  10  CALL GETSPD(MODEL,NPATCH,1,RTMP,ITMP,LTMP,CTMP,IERR)
      IF(IERR.EQ.0) THEN
        GETPCN= RTMP
      ELSEIF(FIRST) THEN
        FIRST= .FALSE.
        IF(NPATCH(4:4).EQ.'*') THEN
          NPATCH(4:4)= '_'
        ELSE
          NPATCH(4:4)= '*'
        ENDIF
        GO TO 10
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C... GXOUTLET deduces the inflow velocity at a fixed pressure boundary
C    from the mass inflow divided by the cell density and area. Allowance
C    is made to relax the inflow velocity.
C    The PATCH commands are:
C    COVAL(patch_name, P1, pco, pval) - standard fixed pressure setting
C    COVAL(patch_name, vel, onlyms, GRND1) - vel is U1 through W2
C    The relaxation is set via SPEDAT as follows:
C      SPEDAT(SET,object_name,RELAX,R,relax) - relaxation factor 0