c
C.... file name GXMISCSO.F 300614
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 'farray'
INCLUDE 'grdloc'
INCLUDE '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==13) THEN
IF(NP=='UP'.AND.(INDVAR==U1.OR.INDVAR==U2)) THEN
CALL FN41(VAL,XU2D,1.0,0.0)
CALL FN25(VAL,-POLRA)
ELSE
IF(NP=='VP'.AND.(INDVAR==V1.OR.INDVAR==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 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
COMMON /UVWCOL/IUC1,IVC1,IWC1,IUVW1(3),IUC2,IVC2,IWC2,IUVW2(23)
COMMON /UVWGBF/ IUC1GCV,IVC1GCV,IWC1GCV,IFILGCV(9)
COMMON/NAMFN/NAMFUN,NAMSUB
LOGICAL FLUID1,MASKPT
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB= 'GXPDRP'
IF(ISC==2 .OR. ISC==3) THEN
L0CO = L0F(CO)
L0VEL= L0F(INDVAR)
CONST= GETPCN(0,IPAT,MASKPT(IPAT))
IF(ISC==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==U1.OR.INDVAR==U2.OR.INDVAR==IUC1GCV.OR.
1 INDVAR==IUC1.OR.INDVAR==IUC2) THEN
LRHON=EAST(LRHO12)
ELSEIF(INDVAR==V1.OR.INDVAR==V2.OR.INDVAR==IVC1GCV.OR.
2 INDVAR==IVC1.OR.INDVAR==IVC2) THEN
LRHON=NORTH(LRHO12)
ELSEIF(INDVAR==W1.OR.INDVAR==W2.OR.INDVAR==IWC1GCV.OR.
2 INDVAR==IWC1.OR.INDVAR==IWC2) 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)>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 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE '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==U1.OR.INDVAR==U2) THEN
LRHON=EAST(LRHO12)
ELSEIF(INDVAR==V1.OR.INDVAR==V2) THEN
LRHON=NORTH(LRHO12)
ELSEIF(INDVAR==W1.OR.INDVAR==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>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 '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==0) THEN
MODEL= 'PDROP_COE'
ELSE
MODEL= 'PDROP_PWR'
ENDIF
10 CALL GETSPD(MODEL,NPATCH,1,RTMP,ITMP,LTMP,CTMP,IERR)
IF(IERR==0) THEN
GETPCN= RTMP
ELSEIF(FIRST) THEN
FIRST= .FALSE.
IF(NPATCH(4:4)=='*') 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 0MAXDMN) THEN
CALL WRITBL
CALL WRITST
CALL WRIT40('Please increase MAXDMN in GXOUTLET')
CALL WRIT2I('Current ',MAXDMN,'; Needed',NUMDMN)
CALL WRITST
CALL SET_ERR(557,'MAXDMN too small in GXINOUT',1)
RETURN
ENDIF
100 CONTINUE
IF(IDMN>1) CALL INDXDM
NOUT(IDMN)=0
NCELL=0
DO IPAT=1,NUMPAT ! loop over patches counting outlets
IF(.NOT.LPTDMN(IPAT)) CYCLE ! not for this domain
IF(IS_OUT(IPAT)) THEN ! count and get no of cells covered
NOUT(IDMN)=NOUT(IDMN)+1
CALL GETPAT(IPAT,IR,GTYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2)
NCELL=MAX(NCELL,(IX2-IX1+1)*(IY2-IY1+1)*(IZ2-IZ1+1))
ENDIF
ENDDO
C... allocate storage
CALL GXMAKE0(L0REL(IDMN),NOUT(IDMN),'RELX') ! relaxation factor
CALL GXMAKE0(L0RAT(IDMN),NOUT(IDMN),'ARAT') ! area ratio
CALL GXMAKE0(L0VEL(IDMN),NOUT(IDMN)*NCELL,'VELIN') ! initial guessed velocity
CALL GXMAKE0(L0OUT(IDMN),NUMPAT,'IOUT') ! sequence number
IF(.NOT.ONEPHS) THEN
CALL GXMAKE0(L0RL2(IDMN),NOUT(IDMN),'RELX2')
CALL GXMAKE0(L0VL2(IDMN),NOUT(IDMN)*NCELL,'VELIN2')
ENDIF
IOUT=0
DO IPAT=1,NUMPAT ! loop again saving settings from Q1
IF(.NOT.LPTDMN(IPAT)) CYCLE
IF(IS_OUT(IPAT)) THEN
COBNAM=OBJNAM(IPAT)
IF(COBNAM.EQ.'NOTSET') COBNAM=NAMPAT(IPAT) ! use patch name if not object
IOUT=IOUT+1
F(L0OUT(IDMN)+IPAT)=IOUT ! save sequence number
VELIN=-999.9
CALL GETSDR(COBNAM,'VELIN',VELIN)
IF(QNE(VELIN,-999.9)) THEN ! set whole patch velocity to initial guess
DO II=1,NCELL
F(L0VEL(IDMN)+(IOUT-1)*NCELL+II)=VELIN
ENDDO
ENDIF
RELX=0.3
CALL GETSDR(COBNAM,'RELAX',RELX) ! relaxation factor
F(L0REL(IDMN)+IOUT)=RELX
IF(.NOT.ONEPHS) THEN
VELIN=-999.9
CALL GETSDR(COBNAM,'VELIN2',VELIN)
IF(QNE(VELIN,-999.9)) THEN
DO II=1,NCELL
F(L0VL2(IDMN)+(IOUT-1)*NCELL+II)=VELIN
ENDDO
ENDIF
RELX=0.3
CALL GETSDR(COBNAM,'RELAX2',RELX)
F(L0RL2(IDMN)+IOUT)=RELX
ENDIF
IF(MASKPT(IPAT)) THEN
CTEMP='^'//NAMPAT(IPAT)
ELSE
CTEMP='!'//NAMPAT(IPAT)
ENDIF
ARAT=1.0
CALL GETSDR('ARATIO',CTEMP,ARAT) ! area ratio
ARAT=AMAX1(ARAT,0.0)
F(L0RAT(IDMN)+IOUT)=ARAT
ENDIF
ENDDO
C... loop back for next FGV domain
IF(IDMN1) THEN
IDMN=1
CALL INDXDM
ENDIF
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C Index for Coefficient - CO
C Index for Value - VAL
ELSEIF(IGR==13.AND.ISC==13) THEN
C------------------- SECTION 14 ------------------- value = GRND1
C For COVAL(patch, velocity, onlyms,GRND1) deduce inflow velocity
C from mass source
IOUT=NINT(F(L0OUT(IDMN)+IPAT)) ! sequence number of outlet
IF(IOUT==0) RETURN ! nothing stored
C... check patch type
CALL GETPAT(IPAT,IREG,RTYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
ITYP=RTYP
C... set area index and flow direction indicator MULT
C MULT = +1 for U at West, V at South, W at Low
C = -1 for U at East, V at North , W at High
C = 0 for all other combinations
IF(ITYP==2) THEN ! East
KGEOM=KXYAE
MULT=ITWO(-1,0,INDVAR==U1.OR.INDVAR==U2)
ELSEIF(ITYP==3) THEN ! West
KGEOM=KXYAW
IF(STORGEO.AND.JXF==1) CALL DOM_BOUND
MULT=ITWO( 1,0,INDVAR==U1.OR.INDVAR==U2)
ELSEIF(ITYP==4) THEN ! North
KGEOM=KXYAN
MULT=ITWO(-1,0,INDVAR==V1.OR.INDVAR==V2)
ELSEIF(ITYP==5) THEN ! South
KGEOM=KXYAS
IF(STORGEO.AND.JYF==1) CALL DOM_BOUND
MULT=ITWO( 1,0,INDVAR==V1.OR.INDVAR==V2)
ELSEIF(ITYP==6) THEN ! High
KGEOM=KXYAH
MULT=ITWO(-1,0,INDVAR==W1.OR.INDVAR==W2)
ELSE ! Low
KGEOM=KXYAL
IF(STORGEO.AND.JZF==1) CALL DOM_BOUND
MULT=ITWO( 1,0,INDVAR==W1.OR.INDVAR==W2)
ENDIF
IF(MULT/=0) THEN ! something to do, velocity normal to face
C... use patch-wise store for mass source, as more reliable
IF(FLUID1(INDVAR)) THEN
KPVMAS=PVMASS ! mass source
L0DEN=L0F(DEN1) ! density
RELX=F(L0REL(IDMN)+IOUT) ! velocity relaxation
L0V=L0VEL(IDMN)+(IOUT-1)*NCELL ! previous sweeps velocity
LBVOUT=LBNAME('VOUT') ! optional storage for inflow velocity
ELSE
KPVMAS=PVMAS2
L0DEN=L0F(DEN2)
RELX=F(L0RL2(IDMN)+IOUT)
L0V=L0VL2(IDMN)+(IOUT-1)*NCELL
LBVOUT=LBNAME('VOU2') ! optional storage for inflow velocity
ENDIF
L0MIN = L0PVAR(KPVMAS,IPAT,IZ)
IF(LBVOUT>0) L0VOUT=L0F(LBVOUT)
C... get area ratio
ARAT=F(L0RAT(IDMN)+IOUT)
C... get index for VAL
L0VAL=L0F(VAL)
RMULT=FLOAT(MULT)
IZADD=(IZ-JZF)*(JXL-JXF+1)*(JYL-JYF+1)
C... High mass sources for W1 are always at IZ+1
IF(LWPZP1) THEN
L0MIN = L0PVAR(KPVMAS,IPAT,IZ+1)
IZADD=(IZ+1-JZF)*(JXL-JXF+1)*(JYL-JYF+1)
IF(LBVOUT>0) L0VOUT=L0F(ANYZ(LBVOUT,IZ+1))
ENDIF
IADD=ITWO(NX*NY,0,INDVAR>6.AND.JZF0) F(L0VOUT+I)=VEL ! store for print/plot
ENDIF
C... set boundary value
F(L0VAL+I)=VEL
ELSE
F(L0VAL+I)=0.0
ENDIF
ENDDO
ENDDO
ELSE
C... nothing to do - velocity is not normal to face
CALL FN1(VAL,0.0)
ENDIF
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
LOGICAL FUNCTION IS_OUT(IPAT)
INCLUDE 'patcmn'
INCLUDE 'objnam'
COMMON /RDATA/ RD1(27),FIXFLU,RD2(32),SAME,RD3(24)
LOGICAL MASKPT,QNE
CHARACTER*14 COBTYP, CTEMP
C... check if object type is outlet
IF(OBJNAM(IPAT)=='NOTSET') THEN
CALL GETCV(IPAT,9,GCO1,GVAL2)
CALL GETCV(IPAT,10,GCO2,GVAL2)
IS_OUT=(GCO1>0.0.AND.QNE(GCO1,FIXFLU)) .OR. (GCO2>0.AND.
1 QNE(GCO2,FIXFLU))
ELSE
IF(MASKPT(IPAT)) THEN
CTEMP='^'//NAMPAT(IPAT)
ELSE
CTEMP='!'//NAMPAT(IPAT)
ENDIF
COBTYP=' '
CALL GETSDC('OBJTYP',CTEMP,COBTYP)
IS_OUT=COBTYP=='OUTLET'.OR.COBTYP=='OPENING'.OR.
1 COBTYP=='APERTURE'
ENDIF
END