c
C.... FILE NAME GXIO.FTN-------------------------------- 260923
C.... FILE NAME GXIO.FTN-------------------------------- 260923
SUBROUTINE GXIO
INCLUDE 'patnos'
INCLUDE 'patcmn'
INCLUDE 'objnam'
INCLUDE 'farray'
INCLUDE 'd_earth/parvar'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
INCLUDE 'd_earth/pbcstr'
INCLUDE 'parear'
INCLUDE 'lbnamer'
COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
COMMON/DRHODP/ITEMP,IDEN/DVMOD/IDVCGR
COMMON/INDAUX/IFL0(4),NSTO,NSOL,IFL(14)
COMMON /INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT,
1 NDOMMX
COMMON /GEODMN0/ I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON /FACES/L0FACE,L0FACZ
COMMON /VRMCMN/L0MSKZ
COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
1 /TSKEMR/ GKTDKP
LOGICAL SLD,QEQ,QNE,FLUID1,DOIT,GRN,LINVRO,MASKPT,QLE,LPTDMN,LEZ,
1 LPAR,EQZ,ISACTIVECELLS
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION STARTS:
C
REAL NORML(3),VELIN(3),VELOUT(3),DIR(3),MASFLOW
PARAMETER (MAXDMN = 10)
CHARACTER*12 COBNAM,COBNM2,BLOCKNAM,CVAL*68,PATNM*8
INTEGER L0PCOEF(MAXDMN),L0BLKOUT(MAXDMN),L0BLKIN(MAXDMN),
1 L0DENIN(MAXDMN),L0TOTVEL(MAXDMN),L0VELI(MAXDMN),
1 L0TOTVELO(MAXDMN),L0VELO(MAXDMN),L0PCOEF2(MAXDMN),
1 L0DENIN2(1),L0TOTVEL2(MAXDMN),L0VELI2(MAXDMN),
1 L0TOTVELO2(MAXDMN),L0VELO2(MAXDMN),L0ARATI(MAXDMN),
1 L0ARATO(MAXDMN),L0KEIN(MAXDMN),L0EPIN(MAXDMN),
1 L0GXO(MAXDMN),L0GXI(MAXDMN),NGXIN(MAXDMN),NGXOUT(MAXDMN),
1 L0TOTAREA(MAXDMN),LINKTO(MAXDMN),LINKFROM(MAXDMN),
1 L0ADDQ(MAXDMN),L0VEL(MAXDMN),L0REL(MAXDMN),L0VL2(MAXDMN),
1 L0RL2(MAXDMN),L0IADQ(MAXDMN),L0TURBIN(MAXDMN),
1 L0IADC(MAXDMN,150),L0ADDC(MAXDMN,150),
1 L0ADMIN(MAXDMN,150),L0ADMAX(MAXDMN,150),L0HUNIT(MAXDMN)
C INTEGER INDPATCH
C REAL GETCOVAL4
SAVE L0PCOEF,L0BLKOUT,L0BLKIN,L0DENIN,L0TOTVEL,L0VELI,
1 L0TOTVELO,L0VELO,L0PCOEF2,L0DENIN2,L0TOTVEL2,L0VELI2,
1 L0TOTVELO2,L0VELO2,L0ARATI,L0ARATO,ISURN,ITM3,
1 L0KEIN,L0EPIN,L0GXO,L0GXI,NGXIN,NGXOUT,L0TOTAREA,LINKTO,
1 LINKFROM,L0ADDQ,L0VEL,L0REL,L0VL2,L0RL2,NCELL,ILIM,L0IADQ,
1 L0TURBIN,L0ADDC,L0IADC,L0ADMIN,L0ADMAX
SAVE LPAR
save l0carea
C
C***********************************************************************
IF(IGR==1.AND.ISC==2) THEN
IF(NUMDMN>MAXDMN) THEN
CALL WRITBL
CALL WRITST
CALL WRIT40('Please increase MAXDMN in GXIO')
CALL WRIT2I('Current ',MAXDMN,'; Needed',NUMDMN)
CALL WRITST
CALL SET_ERR(557,'MAXDMN too small in GXINOUT',1)
RETURN
ENDIF
C.... Loop over all FGV domains
100 CONTINUE
IF(IDMN>1) CALL INDXDM
NGXIN(IDMN)=0
NGXOUT(IDMN)=0
NCELL=0
LPAR=MIMD.AND.NPROC>1
C... Count ins and outs
CALL GXMAKE0(L0GXO(IDMN),NUMPAT,'GXO') ! stores for sequence number
CALL GXMAKE0(L0GXI(IDMN),NUMPAT,'GXI') ! indexed by patch
ILIM=ITWO(GD_NUMPAT,NUMPAT,LPAR) ! for parallel loop over global patches
DO I=1,ILIM
IF(LPAR) THEN
II=GD_INDPAT(I,1); IF(II<0) CYCLE ! if patch not in domain skip
ELSE
II=I
ENDIF
IF(.NOT.LPTDMN(II)) CYCLE ! not this domain
COBNAM=OBJNAM(II); CVAL=' '
CALL GETSDC(COBNAM,'OUTLET',CVAL)
IF(CVAL=='GXO') THEN ! found ANGLED-OUT
NGXOUT(IDMN)=NGXOUT(IDMN)+1
F(L0GXO(IDMN)+II)=NGXOUT(IDMN)
CVAL=' '
CALL GETSDC(COBNAM,'DEDUCED',CVAL) ! external velocity deduced
IF(CVAL=='YES') THEN ! get max number of cells in patch
CALL GETPAT(II,IR,GTYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2)
NCELL=MAX(NCELL,(IX2-IX1+1)*(IY2-IY1+1)*(IZ2-IZ1+1))
ENDIF
CVAL=' '
CALL GETSDC(COBNAM,'DEDUCED2',CVAL) ! external velocity deduced
IF(CVAL=='YES') THEN
CALL GETPAT(II,IR,GTYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2)
NCELL=MAX(NCELL,(IX2-IX1+1)*(IY2-IY1+1)*(IZ2-IZ1+1))
ENDIF
ENDIF
CVAL=' '
CALL GETSDC(COBNAM,'INLET',CVAL)
IF(CVAL=='GXI') THEN ! found ANGLED-IN
NGXIN(IDMN)=NGXIN(IDMN)+1
F(L0GXI(IDMN)+II)=NGXIN(IDMN)
ENDIF
ENDDO
C... Sum number of GXI's over all processors to ensure all have arrays set
IF(LPAR) CALL GLSUMI(NGXIN(IDMN))
IF(NGXOUT(IDMN)>0) THEN ! there are GXOs
NUM=NGXOUT(IDMN) ! stores are indexed by sequence number
CALL GXMAKE0(L0PCOEF(IDMN),NUM,'PCOEFF') ! P1 coefficient
CALL GXMAKE0(L0BLKOUT(IDMN),NUM,'BLKOUT') ! parent blockage
CALL GXMAKE0(L0TOTVELO(IDMN),NUM,'TOTVELO') ! outlet total velocity
CALL GXMAKE0(L0VELO(IDMN),3*NUM,'VELOUT') ! outlet velocity vector
CALL GXMAKE0(L0ARATO(IDMN),NUM,'ARATO') ! outlet area ratio
IF(NCELL>0) THEN
CALL GXMAKE(L0REL(IDMN),NUM,'RELX') ! relaxation for deduced
CALL GXMAKE(L0VEL(IDMN),NUM*NCELL,'VELIN') ! initial guess for deduced
ENDIF
IF(.NOT.ONEPHS) THEN
CALL GXMAKE0(L0PCOEF2(IDMN),NUM,'PCOEFF') ! P2 coefficient
CALL GXMAKE0(L0TOTVELO2(IDMN),NUM,'TOTVELO') ! outlet total velocity
CALL GXMAKE0(L0VELO2(IDMN),3*NUM,'VELOUT') ! outlet velocity vector
IF(NCELL>0) THEN
CALL GXMAKE(L0RL2(IDMN),NUM,'RELX2') ! relaxation for deduced
CALL GXMAKE(L0VL2(IDMN),NUM*NCELL,'VELIN2') ! initial guess for deduced
ENDIF
ENDIF
ENDIF
C
IF(NGXIN(IDMN)>0) THEN ! there are GXIs
NUM=NGXIN(IDMN) ! stores are indexed by sequence number
CALL GXMAKE0(L0BLKIN(IDMN),NUM,'BLKIN') ! parent blockage
CALL GXMAKE0(L0DENIN(IDMN),NUM,'DENIN') ! inlet density
CALL GXMAKE0(L0TOTVEL(IDMN),NUM,'TOTVEL') ! inlet total velocity
CALL GXMAKE0(L0VELI(IDMN),3*NUM,'VELIN') ! inlet velocity vector
CALL GXMAKE0(L0ARATI(IDMN),NUM,'ARATO') ! intlet area ratio
IF(.NOT.ONEPHS) THEN
CALL GXMAKE0(L0DENIN2(IDMN),NUM,'DENIN') ! inlet density
CALL GXMAKE0(L0TOTVEL2(IDMN),NUM,'TOTVEL') ! inlet total velocity
CALL GXMAKE0(L0VELI2(IDMN),3*NUM,'VELIN') ! inlet velocity vector
ENDIF
IF(GRN(ENUT)) THEN
CALL GXMAKE0(L0KEIN(IDMN),NUM,'KEIN')
CALL GXMAKE0(L0EPIN(IDMN),NUM,'EPIN')
CALL GXMAKE0(L0TURBIN(IDMN),NUM,'TURBIN')
ENDIF
CALL GXMAKE0(L0TOTAREA(IDMN),NUM,'TOTAREA')
NUM2=ITWO(ILIM,NUM,LPAR) ! parallel LINKTO is global, sequential is not
CALL GXMAKE0(LINKTO(IDMN),NUM2,'LINKTO')
CALL GXMAKE0(L0ADDQ(IDMN),NUM,'ADDQ')
CALL GXMAKE0(L0IADQ(IDMN),NUM,'IADQ')
DO IPHI=16,NPHI
IF(.NOT.SOLVE(IPHI)) CYCLE
CALL GXMAKE0(L0IADC(IDMN,IPHI),NUM,'IADC')
CALL GXMAKE0(L0ADDC(IDMN,IPHI),NUM,'ADDC')
CALL GXMAKE0(L0ADMIN(IDMN,IPHI),NUM,'MINVAL')
CALL GXMAKE0(L0ADMAX(IDMN,IPHI),NUM,'MAXVAL')
IF(NAME(IPHI)=='MH2O') CALL GXMAKE0(L0HUNIT(IDMN),NUM,
1 'HUNIT')
ENDDO
ENDIF
C
IOUT=0; IIN=0
DO IP=1,ILIM ! loop over patches to get and set data from Q1
IF(LPAR) THEN ! get object name for parallel
COBNAM=GD_OBJNAM(IP)
II=GD_INDPAT(IP,1) ! get local patch number. -ve means not on this proc.
ELSE
II=IP
IF(.NOT.LPTDMN(II)) CYCLE ! not this domain
COBNAM=OBJNAM(II) ! name of object for current patch
ENDIF
CVAL=' '; CALL GETSDC(COBNAM,'OUTLET',CVAL)
IF(II>0.AND.CVAL=='GXO') THEN ! found ANGLED-OUT on this processor
IOUT=IOUT+1
C... deal with angled outlets
PCOEFF=1000. ! default P1 coefficient
CALL GETSDR(COBNAM,'PCOEF',PCOEFF) ! get P1 coefficient
F(L0PCOEF(IDMN)+IOUT)=PCOEFF ! store for later use
BLOCKNAM=' '; IBLK=0
CALL GETSDC(COBNAM,'BLOCKAGE',BLOCKNAM) ! get name of parent
IF(BLOCKNAM/=' ') THEN ! parent name set, so find patch number
DO III=1,ILIM
IF(LPAR) THEN ! get object name for parallel
COBNM2=GD_OBJNAM(III)
ELSE
COBNM2= OBJNAM(III)
ENDIF
IF(COBNM2==BLOCKNAM) THEN
IF(LPAR) THEN
IBLK=GD_INDPAT(III,1)
ELSE
IBLK=III
ENDIF
EXIT
ENDIF
ENDDO
ENDIF
C... If IBLK is zero, angled-out will attach to any blockage it intersects
F(L0BLKOUT(IDMN)+IOUT)=FLOAT(IBLK) ! store for later use
C... Set total outlet velocity if present
TOTVELO=-999.
CALL GETSDR(COBNAM,'TOTVELO',TOTVELO) ! outlet total velocity
F(L0TOTVELO(IDMN)+IOUT)=TOTVELO
C... Set outlet velocity vector if total velocity not set
IF(QEQ(TOTVELO,-999.0)) THEN ! outlet velocity vector
CVAL='NO'
CALL GETSDC(COBNAM,'DEDUCED',CVAL) ! outlet total velocity
IF(CVAL=='NO') THEN ! set external velocity
CALL GETSDR(COBNAM,'UOUT',F(L0VELO(IDMN)+(IOUT-1)*3+1))
! get U
CALL GETSDR(COBNAM,'VOUT',F(L0VELO(IDMN)+(IOUT-1)*3+2))
! get V
CALL GETSDR(COBNAM,'WOUT',F(L0VELO(IDMN)+(IOUT-1)*3+3))
! get W
ELSE ! deduced external velocity
F(L0TOTVELO(IDMN)+IOUT)=-9999.0
VIN=0.0
CALL GETSDR(COBNAM,'VELIN',VIN) ! initial guess
DO IC=1,NCELL
F(L0VEL(IDMN)+(IOUT-1)*NCELL+IC)=VIN
ENDDO
RELX=0.3
CALL GETSDR(COBNAM,'RELAX',RELX) ! relaxation factor
F(L0REL(IDMN)+IOUT)=RELX
ENDIF
ENDIF
C... Set outlet area ratio
ARAT=1.0
CALL GETSDR(COBNAM,'ARATO',ARAT)
F(L0ARATO(IDMN)+IOUT)=ARAT
IF(.NOT.ONEPHS) THEN
PCOEFF=1000000. ! default P2 coefficient
CALL GETSDR(COBNAM,'PCOEF2',PCOEFF) ! get P2 coefficient
F(L0PCOEF2(IDMN)+IOUT)=PCOEFF ! store for later use
C... Set total outlet velocity if present
TOTVELO=-999.
CALL GETSDR(COBNAM,'TOTVELO2',TOTVELO) ! outlet total velocity
F(L0TOTVELO2(IDMN)+IOUT)=TOTVELO
C... Set outlet velocity vector if total velocity not set
IF(QEQ(TOTVELO,-999.0)) THEN ! outlet velocity vector
CVAL='NO'
CALL GETSDC(COBNAM,'DEDUCED2',CVAL) ! outlet total velocity
IF(CVAL=='NO') THEN
CALL GETSDR(COBNAM,'UOUT2',F(L0VELO2(IDMN)+(IOUT-1)*3
1 +1)) ! get U comp
CALL GETSDR(COBNAM,'VOUT2',F(L0VELO2(IDMN)+(IOUT-1)*3
1 +2)) ! get V comp
CALL GETSDR(COBNAM,'WOUT2',F(L0VELO2(IDMN)+(IOUT-1)*3
1 +3)) ! get W comp
ELSE
F(L0TOTVELO2(IDMN)+IOUT)=-9999.0
VIN=0.0
CALL GETSDR(COBNAM,'VELIN2',VIN)
DO IC=1,NCELL
F(L0VL2(IDMN)+(IOUT-1)*NCELL+IC)=VIN
ENDDO
RELX=0.3
CALL GETSDR(COBNAM,'RELAX2',RELX)
F(L0RL2(IDMN)+IOUT)=RELX
ENDIF
ENDIF
ENDIF
C... If STORE(CARE) in Q1, save surface areas into CARE
LBARR=LBNAME('CARE')
IF(LBARR>0) THEN ! need area
C... get patch limits, and sum surface area
L0FACZ0=L0FACZ; IZSTEP=1
IF(II<0) GO TO 10021 ! skip loop for parallel, other processor
CALL GETPAT(II,IREG,TYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
IBLK=NINT(F(L0BLKOUT(IDMN)+IOUT))
DO IZZ=JZF,JZL
IZZM1=(IZZ-1)*NXNY
L0FACZ= L0FACE+IZZM1
L0ARR=L0F(ANYZ(LBARR,IZZ))
IF(MASKPT(II)) L0MSKZ= L0PVAR(22,II,IZZ)
L0PAT=L0PATNO(IDMN)+(IZZ-1)*NXNY ! index for patch-number store
IPW=0
DO IXX=JXF,JXL
DO IYY=JYF,JYL
I=(IXX-1)*NY+IYY; IPW= IPW+ 1
F(L0ARR+I)=0.0
IF(LPAR) THEN
IF(.NOT.ISACTIVECELLS(IXX,IYY,IZZ)) CYCLE
ENDIF
IF(.NOT.SLD(I).AND.LINVRO(IPW))THEN
CALL GET_SURFACE(IXX,IYY,IZZ,IBLK,L0PAT,CAREA,
1 NORML,IPLUS)
F(L0ARR+I)=CAREA
ENDIF
ENDDO
ENDDO
ENDDO
10021 CONTINUE
ENDIF
ENDIF
C
CVAL=' '; CALL GETSDC(COBNAM,'INLET',CVAL)
IF(CVAL=='GXI') THEN ! found ANGLED-IN
C... Deal with Angled inlets, and cartesian velocities in polar
IF(II>0) IIN=IIN+1
C... find parent blockage if set
BLOCKNAM=' '; IBLK=0
CALL GETSDC(COBNAM,'BLOCKAGE',BLOCKNAM) ! get name of parent
IF(BLOCKNAM/=' ') THEN ! parent name set, so find patch number
DO III=1,ILIM
IF(LPAR) THEN ! get object name for parallel
COBNM2=GD_OBJNAM(III)
ELSE
COBNM2= OBJNAM(III)
ENDIF
IF(COBNM2==BLOCKNAM) THEN
IBLK=III
EXIT
ENDIF
ENDDO
ENDIF
IF(II>0) F(L0BLKIN(IDMN)+IIN)=FLOAT(IBLK)
BLOCKNAM=' '; IBLK=0
CALL GETSDC(COBNAM,'LINKTO',BLOCKNAM) ! get name of linked
IF(BLOCKNAM/=' ') THEN
IF(BLOCKNAM=='NEXT') THEN ! find next ANGLED-IN
C... search for next ANGLED-IN in list
DO III=IP+1,ILIM
CVAL=' '
IF(LPAR) THEN
COBNAM=GD_OBJNAM(III)
CALL GETSDC(COBNAM,'INLET',CVAL)
ELSE
CALL GETSDC(OBJNAM(III),'INLET',CVAL)
ENDIF
IF(CVAL=='GXI') THEN
IBLK=III
EXIT
ENDIF
ENDDO
IF(IBLK==0) THEN
CALL WRIT40('No NEXT ANGLED-IN found for '
1 //OBJNAM(II))
ENDIF
ELSEIF(BLOCKNAM=='PREVIOUS') THEN ! find previous ANGLED-IN
C... search for previous ANGLED-IN in list
DO III=IP-1,1,-1
CVAL=' '
IF(LPAR) THEN
COBNAM=GD_OBJNAM(III)
CALL GETSDC(COBNAM,'INLET',CVAL)
ELSE
CALL GETSDC(OBJNAM(III),'INLET',CVAL)
ENDIF
IF(CVAL=='GXI') THEN
IBLK=III
EXIT
ENDIF
ENDDO
IF(IBLK==0) THEN
CALL WRIT40('No PREVIOUS ANGLED-IN found for '
1 //OBJNAM(II))
ENDIF
ELSE ! search for named ANGLED-IN
C... search for named ANGLED-IN
DO III=1,ILIM
IF(LPAR) THEN
COBNAM=GD_OBJNAM(III)
ELSE
COBNAM=OBJNAM(III)
ENDIF
IF(COBNAM==BLOCKNAM) THEN
IBLK=III
EXIT
ENDIF
ENDDO
IF(IBLK==0) THEN
CALL WRIT40('Cannot find linked ANGLED-IN for '
1 //OBJNAM(II))
ENDIF
ENDIF
IF(IBLK==0) THEN ! error finding linked ANGLED-IN
CALL SET_ERR(574,'Cannot find linked ANGLED-IN for '
1 //OBJNAM(II),1)
ENDIF
IF(II>0) THEN
TADD=-999.0
CALL GETSDR(OBJNAM(II),'TADD',TADD) ! get additional temp
IF(QNE(TADD,-999.0)) THEN ! found it
F(L0ADDQ(IDMN)+IIN)=TADD
F(L0IADQ(IDMN)+IIN)=2.
ELSE ! not found, try exit temp
TOUT=-999.0
CALL GETSDR(OBJNAM(II),'TOUT',TOUT) ! get exit temp
IF(QNE(TOUT,-999.0)) THEN ! found it
F(L0ADDQ(IDMN)+IIN)=TOUT
F(L0IADQ(IDMN)+IIN)=1.
ELSE ! not found, try added heat
QADD=0.0
CALL GETSDR(OBJNAM(II),'QADD',QADD) ! get additional heat
F(L0ADDQ(IDMN)+IIN)=QADD
F(L0IADQ(IDMN)+IIN)=0.
ENDIF
ENDIF
IF(ITEM1>0) THEN
TMIN=VARMIN(ITEM1);CALL GETSDR(OBJNAM(II),'TMIN',TMIN)
F(L0ADMIN(IDMN,ITEM1)+IIN)=TMIN
TMAX=VARMAX(ITEM1);CALL GETSDR(OBJNAM(II),'TMAX',TMAX)
F(L0ADMAX(IDMN,ITEM1)+IIN)=TMAX
ELSEIF(SOLVE(H1)) THEN
HMIN=VARMIN(H1); CALL GETSDR(OBJNAM(II),'HMIN',HMIN)
F(L0ADMIN(IDMN,H1)+IIN)=HMIN
HMAX=VARMAX(H1); CALL GETSDR(OBJNAM(II),'HMAX',HMAX)
F(L0ADMAX(IDMN,H1)+IIN)=HMAX
ENDIF
IF(ITEM2>0) THEN
TMIN=VARMIN(ITEM2);CALL GETSDR(OBJNAM(II),'TMIN',TMIN)
F(L0ADMIN(IDMN,ITEM2)+IIN)=TMIN
TMAX=VARMAX(ITEM2);CALL GETSDR(OBJNAM(II),'TMAX',TMAX)
F(L0ADMAX(IDMN,ITEM2)+IIN)=TMAX
ELSEIF(SOLVE(H2)) THEN
HMIN=VARMIN(H2); CALL GETSDR(OBJNAM(II),'HMIN',HMIN)
F(L0ADMIN(IDMN,H2)+IIN)=HMIN
HMAX=VARMAX(H2); CALL GETSDR(OBJNAM(II),'HMAX',HMAX)
F(L0ADMAX(IDMN,H2)+IIN)=HMAX
ENDIF
C.... Source-type flag held in F(L0IADC(IDMN,IPHI),IIN)
C 0 - added flux CFLX
C 1 - fixed exit value COUT
C 2 - added exit value CADD
C 3 - % removal CREM
C IPHI - variable index, II - angled-in counter
DO IPHI=16,NPHI
IF(IPHI==ITEM1.OR.IPHI==ITEM2) CYCLE
IF(.NOT.SOLVE(IPHI)) CYCLE
CADD=-999.0
WRITE(CVAL,'(''CREM'',I3.3)') IPHI; LL=LENGZZ(CVAL)
CALL GETSDR(OBJNAM(II),CVAL(1:LL),CADD) ! Removal %
IF(QNE(CADD,-999.0)) THEN ! found it
F(L0ADDC(IDMN,IPHI)+IIN)=CADD
F(L0IADC(IDMN,IPHI)+IIN)=3.
ELSE ! not found, try additional C
WRITE(CVAL,'(''CADD'',I3.3)') IPHI; LL=LENGZZ(CVAL)
CALL GETSDR(OBJNAM(II),CVAL(1:LL),CADD) ! get additional C
IF(QNE(CADD,-999.0)) THEN ! found it
F(L0ADDC(IDMN,IPHI)+IIN)=CADD
F(L0IADC(IDMN,IPHI)+IIN)=2.
ELSE ! not found, try exit C
COUT=-999.0
WRITE(CVAL,'(''COUT'',I3.3)') IPHI;LL=LENGZZ(CVAL)
CALL GETSDR(OBJNAM(II),CVAL(1:LL),COUT) ! get exit C
IF(QNE(COUT,-999.0)) THEN ! found it
F(L0ADDC(IDMN,IPHI)+IIN)=COUT
F(L0IADC(IDMN,IPHI)+IIN)=1.
ELSE ! not found, try added C flux
CADD=0.0
WRITE(CVAL,'(''CFLX'',I3.3)')IPHI;L=LENGZZ(CVAL)
CALL GETSDR(OBJNAM(II),CVAL(1:L),CADD) ! get additional C
F(L0ADDC(IDMN,IPHI)+IIN)=CADD
F(L0IADC(IDMN,IPHI)+IIN)=0.
ENDIF
ENDIF
ENDIF
WRITE(CVAL,'(''CMIN'',I3.3)') IPHI; LL=LENGZZ(CVAL)
CMIN=VARMIN(IPHI)
CALL GETSDR(OBJNAM(II),CVAL(1:LL),CMIN)
F(L0ADMIN(IDMN,IPHI)+IIN)=CMIN
WRITE(CVAL,'(''CMAX'',I3.3)') IPHI; LL=LENGZZ(CVAL)
CMAX=VARMAX(IPHI)
CALL GETSDR(OBJNAM(II),CVAL(1:LL),CMAX)
F(L0ADMAX(IDMN,IPHI)+IIN)=CMAX
IF(NAME(IPHI)=='MH2O') THEN
IHUNIT=0; CALL GETSDI(OBJNAM(II),'HUNIT',IHUNIT)
F(L0HUNIT(IDMN)+IIN) = IHUNIT
ENDIF
ENDDO
ENDIF
ENDIF
IINL=ITWO(IP,IIN,LPAR)
F(LINKTO(IDMN)+IINL)=FLOAT(IBLK)
C... set inlet density
IF(RHO1<=0.0) THEN ! initialise inlet density
RHOIN=1.0
ELSE
RHOIN=RHO1
ENDIF
CALL GETSDR(COBNAM,'DENIN',RHOIN) ! get density at inlet
IF(II>0) F(L0DENIN(IDMN)+IIN)=RHOIN
C... Set inlet area ratio
ARAT=1.0
CALL GETSDR(COBNAM,'ARATI',ARAT)
IF(II>0) F(L0ARATI(IDMN)+IIN)=ARAT
C... Set total inlet velocity if present
TOTVEL=-999.
CALL GETSDR(COBNAM,'TOTVEL',TOTVEL) ! inlet total velocity
IF(QNE(TOTVEL,-999.0).and.ii>0) THEN
C... for -ve total velocity, set density to GRND, to flag in-cell density
IF(TOTVEL<0.0) F(L0DENIN(IDMN)+IIN)=GRND
ENDIF
IF(II>0) F(L0TOTVEL(IDMN)+IIN)=TOTVEL
C... Set inlet velocity vector if total velocity not set
IF(QEQ(TOTVEL,-999.0)) THEN ! inlet velocity vector
UVEL=0.0; VVEL=0.0; WVEL=0.0
CALL GETSDR(COBNAM,'UIN',UVEL) ! get U cmp
CALL GETSDR(COBNAM,'VIN',VVEL) ! get V cmp
CALL GETSDR(COBNAM,'WIN',WVEL) ! get W cmp
IF(II>0) THEN ! on this processor
F(L0VELI(IDMN)+(IIN-1)*3+1)=UVEL
F(L0VELI(IDMN)+(IIN-1)*3+2)=VVEL
F(L0VELI(IDMN)+(IIN-1)*3+3)=WVEL
ENDIF
ENDIF
C... Set inlet volumetric or mass flow rate
CALL SUB2R(VOLFLOW,-999.0,MASFLOW,-999.0)
CALL GETSDR(COBNAM,'VOLFLOW',VOLFLOW) ! inlet volume flow
CALL GETSDR(COBNAM,'MASFLOW',MASFLOW) ! inlet mass flow
C... If STORE(CARE) in Q1, save surface ares into CARE
LBARR=LBNAME('CARE')
IF(QNE(VOLFLOW,-999.0).OR.QNE(MASFLOW,-999.0).OR.
1 GRN(ENUT).OR.ISG62==1.OR.LBARR>0) THEN ! need area
C... get patch limits, and sum surface area
TOTAREA=0.0
L0FACZ0=L0FACZ
IF(II<0) GO TO 1002 ! skip loop for parallel, other processor
CALL GETPAT(II,IREG,TYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
IBLK=NINT(F(L0BLKIN(IDMN)+IIN))
DO IZZ=JZF,JZL
IZZM1=(IZZ-1)*NXNY
L0FACZ= L0FACE+IZZM1
IF(LBARR>0) L0ARR=L0F(ANYZ(LBARR,IZZ))
IF(MASKPT(II)) L0MSKZ= L0PVAR(22,II,IZZ)
L0PAT=L0PATNO(IDMN)+(IZZ-1)*NXNY ! index for patch-number store
IPW=0
DO IXX=JXF,JXL
DO IYY=JYF,JYL
I=(IXX-1)*NY+IYY; IPW= IPW+ 1
IF(LBARR>0) F(L0ARR+I)=0.0
IF(LPAR) THEN
IF(.NOT.ISACTIVECELLS(IXX,IYY,IZZ)) CYCLE
ENDIF
IF(.NOT.SLD(I).AND.LINVRO(IPW))THEN
CALL GET_SURFACE(IXX,IYY,IZZ,IBLK,L0PAT,CAREA,
1 NORML,IPLUS)
TOTAREA=TOTAREA+CAREA
IF(LBARR>0) F(L0ARR+I)=CAREA
ENDIF
ENDDO
ENDDO
ENDDO
1002 CONTINUE
IF(LPAR) CALL GLSUM(TOTAREA) ! sum over all processors
IF(II>0) F(L0TOTAREA(IDMN)+IIN)=TOTAREA
C.... deduce total velocity from flowrate/area
IF(QNE(VOLFLOW,-999.0).AND.II>0) THEN
F(L0TOTVEL(IDMN)+IIN)=VOLFLOW/(ARAT*TOTAREA+TINY)
C... for -ve volume flow rate, set density to GRND, to flag in-cell density
IF(VOLFLOW<0.0) F(L0DENIN(IDMN)+IIN)=GRND
ELSEIF(QNE(MASFLOW,-999.0).AND.II>0) THEN
F(L0TOTVEL(IDMN)+IIN)=MASFLOW/(ARAT*RHOIN*TOTAREA+TINY)
C... for mass flow rate, set density to -ve, to flag mass flow
F(L0DENIN(IDMN)+IIN)=-RHOIN
ENDIF
L0FACZ=L0FACZ0
ENDIF
IF(GRN(ENUT).AND.II>0) THEN
C... Inlet values for Turbulence
TURBIN=-999.; CALL GETSDR(COBNAM,'TURBIN',TURBIN)
F(L0TURBIN(IDMN)+IIN)=TURBIN
IF(QNE(TURBIN,-999.0)) THEN ! Intensity was set
IF(QEQ(F(L0TOTVEL(IDMN)+IIN),-999.)) THEN
VEL= (F(L0VELI(IDMN)+(IIN-1)*3+1)**2+
1 F(L0VELI(IDMN)+(IIN-1)*3+2)**2+
1 F(L0VELI(IDMN)+(IIN-1)*3+3)**2)**0.5
ELSE
VEL=F(L0TOTVEL(IDMN)+IIN)
ENDIF
GKEIN=(TURBIN*VEL)**2
F(L0KEIN(IDMN)+IIN)=GKEIN
RHYDR=(TOTAREA/(4.*ATAN(1.)))**0.5 ! hydraulic radius
F(L0EPIN(IDMN)+IIN)=CMUCD**0.75*GKEIN**1.5/(0.1*RHYDR)
ELSE ! Intensity not set, so get inlet KE and EP
GKEIN=-999.; CALL GETSDR(COBNAM,'KEIN',GKEIN)
IF(QNE(GKEIN,-999.)) THEN
IF(SOLVE(KE)) THEN
F(L0KEIN(IDMN)+IIN)=GKEIN
ENDIF
ENDIF
GEPIN=-999.; CALL GETSDR(COBNAM,'EPIN',GEPIN)
IF(QNE(GEPIN,-999.)) THEN
IF(SOLVE(EP).OR.SOLVE(LBOMEG).OR.SOLVE(LBET).OR.
1 SOLVE(LBKP) .OR.SOLVE(LBKT)) THEN
F(L0EPIN(IDMN)+IIN)=GEPIN
ENDIF
ENDIF
ENDIF
ENDIF
IF(.NOT.ONEPHS) THEN
C... set inlet density
IF(RHO2<=0.0) THEN ! initialise inlet density
RHOIN=1000.0
ELSE
RHOIN=RHO2
ENDIF
CALL GETSDR(COBNAM,'DENIN2',RHOIN) ! get density at inlet
IF(II>0) F(L0DENIN2(IDMN)+IIN)=RHOIN
C... Set total inlet velocity if present
TOTVEL=-999.
CALL GETSDR(COBNAM,'TOTVEL2',TOTVEL) ! inlet total velocity
IF(QNE(TOTVEL,-999.0).AND.II>0) THEN
C... for -ve total velocity, set density to -1, to flag in-cell density
IF(TOTVEL<0.0) F(L0DENIN2(IDMN)+IIN)=-1.0
ENDIF
IF(II>0) F(L0TOTVEL2(IDMN)+IIN)=TOTVEL
C... Set inlet velocity vector if total velocity not set
IF(QEQ(TOTVEL,-999.0)) THEN ! inlet velocity vector
UVEL=0.0; VVEL=0.0; WVEL=0.0
CALL GETSDR(COBNAM,'UIN2',UVEL) ! get U
CALL GETSDR(COBNAM,'VIN2',VVEL) ! get V
CALL GETSDR(COBNAM,'WIN2',WVEL) ! get W
IF(II>0) THEN
F(L0VELI2(IDMN)+(IIN-1)*3+1)=UVEL
F(L0VELI2(IDMN)+(IIN-1)*3+2)=VVEL
F(L0VELI2(IDMN)+(IIN-1)*3+3)=WVEL
ENDIF
ENDIF
C... Set inlet volumetric flow rate
CALL SUB2R(VOLFLOW,-999.0,MASFLOW,-999.0)
CALL GETSDR(COBNAM,'VOLFLOW2',VOLFLOW) ! inlet volume flow
CALL GETSDR(COBNAM,'MASFLOW2',MASFLOW) ! inlet mass flow
IF(QNE(VOLFLOW,-999.0).OR.QNE(MASFLOW,-999.0)) THEN ! extraction rate set
C... get patch limits, and sum surface area
TOTAREA=0.0
IF(II<0) GO TO 1003
CALL GETPAT(II,IREG,TYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,
1 JTL)
L0FACZ0=L0FACZ
DO IZZ=JZF,JZL
IZZM1=(IZZ-1)*NXNY
L0FACZ= L0FACE+IZZM1
IF(MASKPT(II)) L0MSKZ= L0PVAR(22,II,IZZ)
L0PAT=L0PATNO(IDMN)+(IZZ-1)*NXNY ! index for patch-number store
IPW=0
IBLK=NINT(F(L0BLKIN(IDMN)+IIN))
DO IXX=JXF,JXL
DO IYY=JYF,JYL
I=(IXX-1)*NY+IYY
IPW= IPW+ 1
IF(.NOT.SLD(I).AND.LINVRO(IPW)) THEN
CALL GET_SURFACE(IXX,IYY,IZZ,IBLK,L0PAT,
1 CAREA,NORML,IPLUS)
TOTAREA=TOTAREA+CAREA
ENDIF
ENDDO
ENDDO
ENDDO
1003 CONTINUE
IF(LPAR) CALL GLSUM(TOTAREA)
IF(QNE(VOLFLOW,-999.0).AND.II>0) THEN
C.... deduce total velocity from flowrate/area
F(L0TOTVEL2(IDMN)+IIN)=VOLFLOW/(ARAT*TOTAREA+TINY)
C... for -ve flow rate, set density to -1, to flag in-cell density
IF(VOLFLOW<0.0) F(L0DENIN2(IDMN)+IIN)=-1.0
ELSEIF(QNE(MASFLOW,-999.0).AND.II>0) THEN
F(L0TOTVEL2(IDMN)+IIN)=MASFLOW/(ARAT*TOTAREA*RHOIN+
1 TINY)
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
C... loop back for next FGV domain
IF(IDMN1) THEN
IDMN=1
CALL INDXDM
ENDIF
C... set index for SURN
ISURN=LBSURN
C... set index for T3
ITM3=LBNAME('T3')
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==3) THEN
C------------------- SECTION 3 ------------- coefficient = GRND2
C
C Fixed pressure at OPENING - set pressure
C ----------------------------------------
C The patch name can be anything. The COefficient for P1 must be set to GRND2.
C The COefficient for T3 can also be set to GRND2 to set a radiative loss
C to the surroundings. The VALue for T3 sets the external temperature.
C
C Inputs:
C SPEDAT(SET,object_name,OUTLET,C,GXO) - identifies this as an Angled-out patch
C SPEDAT(SET,object_name,BLOCKAGE,C,name_of_blockage) - sets the
C name of the blockage object this outlet sits on. If not set,
C it will sit on any blockage within the bounding box.
C SPEDAT(SET,object_name,PCOEF,R,pressure_coefficient) - sets
C the coefficient used to fix the pressure. If not set, 1000 is used.
C
IOUT=NINT(F(L0GXO(IDMN)+IPAT))
IIN=NINT(F(L0GXI(IDMN)+IPAT))
DOIT=(IOUT+IIN)>0
IF(DOIT) THEN
IPT=MAX(IOUT,IIN)
DOIT=.FALSE.
IF(INDVAR==P1.OR.INDVAR==P2.OR.INDVAR==R1.OR.
1 INDVAR==R2) THEN
C... Set CO for fixed pressure
IF(FLUID1(INDVAR)) THEN
PCOEFF=F(L0PCOEF(IDMN) +IPT) ! get P1 coefficient
ELSE
PCOEFF=F(L0PCOEF2(IDMN)+IPT) ! get P2 coefficient
ENDIF
IF(LEZ(PCOEFF)) THEN
CALL GETCOV(NPATCH,INDVAR,GCO,GVAL)
L0P1=L0F(P1)
ENDIF
DOIT=.TRUE.
ELSEIF(INDVAR==ITM3) THEN
C... prepare CO for T3
SIGMA=5.6697E-08
L0T3=L0F(ITM3)
CALL GETCOV(NPATCH,INDVAR,GCO,GVAL)
TW =GVAL+TEMP0
TW2=TW**2
TW3=TW**3
DOIT=.TRUE.
ENDIF
IF(DOIT) THEN
L0CO=L0F(CO) ! index for COefficient
L0PAT=L0PATNO(IDMN)+(IZ-1)*NXNY ! index for patch-number store
IF(HOL.OR.SURF.OR.LEZ(PCOEFF)) THEN
IDEN=ITWO(DEN1,DEN2,FLUID1(INDVAR))
L0DEN=L0F(IDEN)
ENDIF
C... get parent blockage patch number
IF(IOUT>0) THEN
IBLK=NINT(F(L0BLKOUT(IDMN)+IPT))
ELSE
IBLK=NINT(F(L0BLKIN(IDMN)+IPT))
ENDIF
IPW=0
C... loop over all cells in bounding box
DO IX=IXF,IXL
DO IY=IYF,IYL
I=(IX-1)*NY+IY
F(L0CO+I)=0.0 ! initialise CO to 0.0
C... locate cells which are in fluid, but have patch number >0,
C i.e. lie on fluid-side of intersection of outlet and blockage
IPW=IPW+1
IF(.NOT.SLD(I).AND.LINVRO(IPW)) THEN
CALL GET_SURFACE(IX,IY,IZ,IBLK,L0PAT,CAREA,NORML,IPLS)
C... CAREA is surface area in affected cell.
IF(CAREA>0.0) THEN
IF(INDVAR==ITM3) THEN
C... Set CO for T3 - equivalent of *RAD patch
T3=F(L0T3+I)+TEMP0
F(L0CO+I)=SIGMA*CAREA*(T3**3+T3**2*TW+T3*TW2+TW3)
ELSE
C... Set CO for mass flow - Pcoeff*area
IF(HOL.OR.SURF) THEN
C... For HOL or SURF, CO is set to local density*surface_area_in_cell
F(L0CO+I)=F(L0DEN+I)*CAREA
ELSE
IF(LEZ(PCOEFF)) THEN
C... Mass flow is m = A.sqrt(abs(C)(V-Pp))
GCO=ABS(PCOEFF)*F(L0DEN+I)
GCO=GCO/(SQRT(ABS(GCO*(GVAL-F(L0P1+I))))+TINY)
F(L0CO+I)=GCO*CAREA
ELSE
C... otherwise CO is set to PCOEF*surface_area_in_cell
F(L0CO+I)=PCOEFF*CAREA
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDIF
ENDIF
ELSEIF(IGR==13.AND.ISC==14) THEN
C------------------- SECTION 14 ------------------- value = GRND2
C
C Fixed pressure at OPENING - set external velocities
C -------------------------
C
C The patch name can be anything. The VALue for velocity must be set to GRND2
C if it is to be set to a user-set or deduced value.
C
C Inputs:
C SPEDAT(SET,object_name,OUTLET,C,GXO) - identifies this as an Angled-out patch
C SPEDAT(SET,object_name,TOTVELO,R,velocity_normal_to_surface) - this
C sets the velocity magnitude normal to the blockage-object surface.
C If it is not set, the following SPEDATs will be used:
C SPEDAT(SET,object_name,UOUT,R,U_component)
C SPEDAT(SET,object_name,VOUT,R,V_component)
C SPEDAT(SET,object_name,WOUT,R,W_component)
C These set the components of the vector at the outlet.
C Alternatively, the velocity can be deduced from mass-flux/density.
C SPEDAT(SET,object_name,DEDUCED,C,YES)
C SPEDAT(SET,object_name,VELIN,R,initial_guess_for_external_velocity)
C SPEDAT(SET,object_name,RELAX,R,linear_relaxation_factor_for_external_velocity)
C
IPT=NINT(F(L0GXO(IDMN)+IPAT))
IF(IPT>0) THEN ! current patch is for ANGLED-OUT
L0PAT=L0PATNO(IDMN)+(IZ-1)*NXNY ! index for patch-number store
L0VAL=L0F(VAL) ! index for VALue
IF(FLUID1(INDVAR)) THEN
TOTVELO=F(L0TOTVELO(IDMN) +IPT) ! total velocity at outlet
L0VOUT=L0VELO(IDMN)
KPVMAS=PVMASS
L0DEN=L0F(DEN1)
IF(TOTVELO==-9999.0.AND.NCELL>0) THEN ! deduced
RELX=F(L0REL(IDMN)+IPT)
L0V=L0VEL(IDMN)+(IPT-1)*NCELL
ENDIF
LBVOUT=LBNAME('VOUT')
ELSE
TOTVELO=F(L0TOTVELO2(IDMN)+IPT) ! total velocity at outlet
L0VOUT=L0VELO2(IDMN)
KPVMAS=PVMAS2
L0DEN=L0F(DEN2)
IF(TOTVELO==-9999.0.AND.NCELL>0) THEN ! deduced
RELX=F(L0RL2(IDMN)+IPT)
L0V=L0VL2(IDMN)+(IPT-1)*NCELL
ENDIF
LBVOUT=LBNAME('VOU2')
ENDIF
IF(LBVOUT>0) L0VOUT=L0F(LBVOUT)
ARAT=F(L0ARATO(IDMN)+IPT) ! area ratio
IF(QEQ(TOTVELO,-999.)) THEN ! velocity components set
VELOUT(1)=F(L0VOUT+(IPT-1)*3+1) ! get U component at outlet
VELOUT(2)=F(L0VOUT+(IPT-1)*3+2) ! get V component at outlet
VELOUT(3)=F(L0VOUT+(IPT-1)*3+3) ! get W component at outlet
IF(INDVAR==U1.OR.INDVAR==U2) THEN ! set direction vector
CALL VECTOR(DIR,1.,0.,0.); IDIR=1
ELSEIF(INDVAR==V1.OR.INDVAR==V2) THEN
CALL VECTOR(DIR,0.,1.,0.); IDIR=2
ELSEIF(INDVAR==W1.OR.INDVAR==W2) THEN
CALL VECTOR(DIR,0.,0.,1.); IDIR=3
ELSE
IDIR=0
ENDIF
IF(INDVAR>=U1.AND.INDVAR<=W2) VEL=DOT(VELOUT,DIR) ! get velocity magnitude in velocity d
ELSEIF(QEQ(TOTVELO,-9999.).OR.QEQ(TOTVELO,SAME)) THEN ! velocity is deduced
L0MIN = L0PVAR(KPVMAS,IPAT,IZ)
L0MINH= L0PVAR(KPVMAS,IPAT,IZ+1)
IDIR=ITWO((INDVAR-1)/2,0,INDVAR>=U1.AND.INDVAR<=W2) ! for velocity IDIR is based on dire
ENDIF
L0XU=L0F(XU2D)
L0XG=L0F(XG2D)
IBLK=NINT(F(L0BLKOUT(IDMN)+IPT)) ! get parent blockage patch number
IPW=0
IF(TOTVELO==-9999.0.AND.NCELL>0) THEN ! deduced
CALL GETPAT(IPAT,IREG,RTYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
IZADD=(IZ-JZF)*(JXL-JXF+1)*(JYL-JYF+1)
ENDIF
C... Deduce external value from VFOL
IF(INDVAR==ISURN) THEN ! Phase 1
L0SURN=L0F(ISURN)
L0MIN = L0PVAR(KPVMAS,IPAT,IZ)
CALL GETCV(IPAT,IVFOL,GCO,GVAL)
RHOEXT=1./GVAL ! VAL for VFOL has been set to 1/density
RHOLIQ=F(INDPRTB(IPRPSA,0)+1)
SURNEXT=RTWO(1.,0., ABS((RHOEXT-RHOLIQ)/RHOLIQ)<=1E-3)
ELSEIF(INDVAR==LBSRN2) THEN ! Phase 3
L0SURN=L0F(LBSRN2)
L0MIN = L0PVAR(KPVMAS,IPAT,IZ)
CALL GETCV(IPAT,LBVFL2,GCO,GVAL)
RHOEXT=1./GVAL ! VAL for VFOL has been set to 1/density
RHOLIQ=F(INDPRTB(IPRPSC,0)+1)
SURNEXT=RTWO(1.,0., ABS((RHOEXT-RHOLIQ)/RHOLIQ)<=1E-3)
ELSEIF(SURF) THEN ! not SURN / SRN2 but VOF
CALL GETCV(IPAT,IVFOL,GCO,GVAL)
RHOEXT=1./(GVAL+TINY) ! VAL for VFOL has been set to 1/density
RHOLIQ=F(INDPRTB(IPRPSA,0)+1)
SURNEXT=RTWO(1.,0., ABS((RHOEXT-RHOLIQ)/RHOLIQ)<=1E-3)
IF(LBVFL2>0) THEN ! Phase 3
CALL GETCV(IPAT,LBVFL2,GCO,GVAL)
RHOXT2=1./(GVAL+TINY) ! VAL for VFOL has been set to 1/density
RHOLQ2=F(INDPRTB(IPRPSC,0)+1)
SURNXT2=RTWO(1.,0., ABS((RHOXT2-RHOLQ2)/RHOLQ2)<=1E-3)
ELSE ! no Phase 3
SURNXT2=0.0
ENDIF
IF(SURNEXT==0.0.AND.SURNXT2==0.0) THEN ! No Phase 1 or 3, must be 2
RHOEXT=F(INDPRTB(IPRPSB,0)+1)
ELSEIF(SURNEXT==1.0) THEN ! Phase 1
RHOEXT=F(INDPRTB(IPRPSA,0)+1)
ELSE ! Phase 3
RHOEXT=F(INDPRTB(IPRPSC,0)+1)
ENDIF
ENDIF
DO IX=IXF,IXL ! loop over all cells in bounding box
DO IY=IYF,IYL
I=(IX-1)*NY+IY
IPW=IPW+1
IJKDM=(IX-1)*NY+IY+(IZ-1)*NXNY ! cell index in 3D
F(L0VAL+I)=0.0 ! initialise VAL to 0.0
C... locate cells which are in fluid, but have patch number >0,
C i.e. lie on fluid-side of intersection of inlet and blockage
IF(.NOT.SLD(I).AND.LINVRO(IPW)) THEN
CALL GET_SURFACE(IX,IY,IZ,IBLK,L0PAT,CAREA,NORML,IPLUS)
C... CAREA is surface area in affected cell.
IF(CAREA>0.0) THEN
C... VAL set to velocity component in relevant direction
IF(QEQ(TOTVELO,-999.)) THEN ! velocity vector has been set
IF(CARTES.OR.INDVAR>V2) THEN
IF(INDVAR==ISURN.OR.INDVAR==LBSRN2) THEN
VELIO=ABS(DOT(VELOUT,NORML))
ELSE
VELIO=VEL
ENDIF
ELSE ! Polar U1-V2
UC=VELOUT(1) ! set cartesian components
VC=VELOUT(2)
IF(INDVAR=U1.AND.INDVAR<=W2.AND.ISWEEP>FSWEEP+1)
1 THEN
L0MASS=ITWO(L0MIN,L0MINH,IPLUS==0)
IF(SURF) THEN ! For vof use density of incoming fluid
VEL=F(L0MASS+IPW)/(ARAT*CAREA*RHOEXT+TINY)
ELSE ! otherwise use in-cell density
VEL=F(L0MASS+IPW)/(ARAT*CAREA*F(L0DEN+I)+TINY)
ENDIF
IF(NCELL>0) THEN
IV=IPW+IZADD
VEL=(1.-RELX)*F(L0V+IV)+RELX*VEL ! relax inflow velocity
F(L0V+IV)=VEL
ENDIF
IF(LBVOUT>0) F(L0VOUT+I)=VEL
ELSE ! not velocity, used stored outflow velocity
IF(NCELL>0) THEN
IV=IPW+IZADD
VEL=F(L0V+IV)
ENDIF
ENDIF
IF(IDIR>0) THEN ! for velocity, resolve into surface direction
VELIO=NORML(IDIR)*VEL
ELSE
VELIO=VEL
ENDIF
ELSEIF(QEQ(TOTVELO,SAME)) THEN
FLO=F(L0MIN+IPW) ! mass flow in patch, +ve in, -ve out
IF(FLO>=0.0) THEN ! inflow
VELIO=FLO/(F(L0DEN+I)*CAREA*ARAT+TINY)
ENDIF
ELSE ! velocity normal to surface has been set
IF(INDVAR>=U1.AND.INDVAR<=W2) THEN
IDIR=(INDVAR-1)/2
VELIO=TOTVELO*NORML(IDIR) ! resolve into grid direction
ELSE
VELIO=TOTVELO
ENDIF
ENDIF
IF(INDVAR>=U1.AND.INDVAR<=W1) THEN ! momentum sources
F(L0VAL+I)=VELIO
ELSEIF(INDVAR==ISURN.OR.INDVAR==LBSRN2) THEN
FLO=F(L0MIN+IPW) ! mass flow in patch, +ve in, -ve out
IF(FLO<=0.0) THEN ! outflow
F(L0VAL+I)=FLO*F(L0SURN+I)/(F(L0DEN+I)) ! take in-cell value
ELSE ! inflow
F(L0VAL+I)=SURNEXT*ABS(VELIO)*CAREA ! take external value
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDIF
C
C Fixed flow rate at Angled INLET
C -------------------------------
C
C The patch name can be anything. The COVALs should be:
C COVAL(GXIn, P1, FIXFLU, GRND2)
C COVAL(GXIn, vel, ONLYMS, GRND2)
C
C Inputs:
C SPEDAT(SET,object_name,INLET,C,GXI) - identifies this as an Angled-in patch
C SPEDAT(SET,object_name,BLOCKAGE,C,name_of_blockage) - sets the
C name of the blockage object this inlet sits on. If not set,
C it will sit on any blockage within the bounding box.
C SPEDAT(SET,object_name,TOTVEL,R,velocity_normal_to_surface) - this
C sets the velocity magnitude normal to the blockage-object surface.
C If it is not set, the following SPEDATs will be used:
C SPEDAT(SET,object_name,UIN,R,U_component)
C SPEDAT(SET,object_name,VIN,R,V_component)
C SPEDAT(SET,object_name,WIN,R,W_component)
C These set the components of the inlet vector.
C or SPEDAT(SET,object_name,VOLFLOW,R,volumetric_flow_rate) - this sets
C the volumetric flow rate. The velocity is deduced from flowrate
C divided by the total inlet area.
C or SPEDAT(SET,object_name,MASFLOW,R,mass_flow_rate) - this sets
C the mass flow rate. The velocity is deduced from flowrate
C divided by the total inlet area divided by the inlet density.
C SPEDAT(SET,object_name,DENIN,R,inlet_density) - this sets the inlet
C density. if not set, RHO1 is used. If that is GRNDn, 1.0 is used.
C For fixed outflows, the inlet density is stored as -1, to flag the
C use of the in-cell density.
C
C... To link the outflow from one ANGELD-IN to the inflow to another, use
C SPEDAT(SET,object_name, LINKTO,C,previous / next)
C previous links to the immediately preceeding ANGLED-IN,
C next links to the following ANGLED-IN
IPT=NINT(F(L0GXI(IDMN)+IPAT)) ! angled-in sequence number
IF(IPT>0) THEN ! current patch is for ANGLED-IN
L0PAT=L0PATNO(IDMN)+(IZ-1)*NXNY ! index for patch-number store
L0VAL=L0F(VAL) ! index for VALue
IF(FLUID1(INDVAR)) THEN
RHOIN=F(L0DENIN(IDMN) +IPT) ! get density at inlet
C... RHOin=GRND signals -ve normal velocity or -ve volume flow rate
C... RHOin < 0 signals fixed mass flow rate
TOTVEL=F(L0TOTVEL(IDMN)+IPT) ! total velocity at inlet
L0VIN=L0VELI(IDMN)
L0DEN=L0F(DEN1) ! index for density
ELSE
RHOIN=F(L0DENIN2(IDMN) +IPT) ! get density at inlet
TOTVEL=F(L0TOTVEL2(IDMN)+IPT) ! total velocity at inlet
L0VIN=L0VELI2(IDMN)
L0DEN=L0F(DEN2) ! index for density
ENDIF
ARAT=F(L0ARATI(IDMN)+IPT) ! area ratio
IF(QEQ(TOTVEL,-999.)) THEN
VELIN(1)=F(L0VIN+(IPT-1)*3+1) ! get U component at inlet
VELIN(2)=F(L0VIN+(IPT-1)*3+2) ! get V component at inlet
VELIN(3)=F(L0VIN+(IPT-1)*3+3) ! get W component at inlet
IF(INDVAR==U1.OR.INDVAR==U2) THEN ! set direction vector
CALL VECTOR(DIR,1.,0.,0.)
ELSEIF(INDVAR==V1.OR.INDVAR==V2) THEN
CALL VECTOR(DIR,0.,1.,0.)
ELSEIF(INDVAR==W1.OR.INDVAR==W2) THEN
CALL VECTOR(DIR,0.,0.,1.)
ENDIF
VEL=DOT(VELIN,DIR) ! get velocity magnitude in velocity direction
ENDIF
IBLK=NINT(F(L0BLKIN(IDMN)+IPT))
IDIR=(INDVAR-1)/2
IPW=0
l0vall=0; l0inv=0
if(indvar==u1) then
if(lbname('UVAL')>0) then
l0vall=l0f(lbname('UVAL'))
endif
elseif(indvar==v1) then
if(lbname('VVAL')>0) then
l0vall=l0f(lbname('VVAL'))
endif
elseif(indvar==w1) then
if(lbname('WVAL')>0) then
l0vall=l0f(lbname('WVAL'))
endif
elseif(indvar==p1.or.indvar==r1) then
if(lbname('PVAL')>0) then
l0vall=l0f(lbname('PVAL'))
endif
if(lbname('#INV')>0) then
l0inv=l0f(lbname('#INV'))
endif
elseif(indvar==u2) then
if(lbname('UVL2')>0) then
l0vall=l0f(lbname('UVL2'))
endif
elseif(indvar==v2) then
if(lbname('VVL2')>0) then
l0vall=l0f(lbname('VVL2'))
endif
elseif(indvar==w2) then
if(lbname('WVL2')>0) then
l0vall=l0f(lbname('WVL2'))
endif
elseif(indvar==p2.or.indvar==r2) then
if(lbname('PVL2')>0) then
l0vall=l0f(lbname('PVL2'))
endif
endif
IF(.NOT.CARTES) THEN
L0XU=L0F(XU2D)
L0XG=L0F(XG2D)
IF(IURVAL/=0) L0RAD=L0F(RG2D)
ENDIF
IF(INDVAR==ISURN.OR.INDVAR==LBSRN2) THEN
ISRN=ITWO(ISURN,LBSRN2,INDVAR==ISURN)
IVFL=ITWO(IVFOL,LBVFL2,INDVAR==ISURN)
CALL GETCV(IPAT,IVFL,GCO,GVFOL)
RHOEXT=1./GVFOL ! VAL for VFOL has been set to 1/density
IPRPSAC=ITWO(IPRPSA,IPRPSC,INDVAR==ISURN)
RHOLIQ=F(INDPRTB(IPRPSAC,0)+1) ! get fluid density, heavy or third
SURNEXT=RTWO(1.,0., ABS((RHOEXT-RHOLIQ)/RHOLIQ)<=1E-3) ! 1 if fluid, 0 if not
ENDIF
DO IX=IXF,IXL ! loop over all cells in bounding box
DO IY=IYF,IYL
I=(IX-1)*NY+IY
IPW=IPW+1
F(L0VAL+I)=0.0 ! initialise VAL to 0.0
C... locate cells which are in fluid, but have patch number >0,
C i.e. lie on fluid-side of intersection of inlet and blockage
CAREA=0.0 ! initialise
if(l0inv>0) then
f(l0inv+i)=rtwo(1.,0.,linvro(ipw))
endif
IF(.NOT.SLD(I).AND.LINVRO(IPW)) THEN
C... In fluid, and within angled-in object. Check if any surface in cell
CALL GET_SURFACE(IX,IY,IZ,IBLK,L0PAT,CAREA,NORML,IPLUS)
ELSEIF(INDVAR==W1.OR.INDVAR==W2) THEN
C... Didn't find any surface. For W, look 1 slab up in case area there
L0MSKZ0=L0MSKZ ! save
IF(MASKPT(IPAT))L0MSKZ=L0PVAR(22,IPAT,IZ+1) ! advance 1 slab
IF(.NOT.SLD(I+NXNY).AND.LINVRO(IPW)) THEN
C... Check if surface on next slab
CALL GET_SURFACE(IX,IY,IZ+1,IBLK,L0PAT,CAREA,NORML,
1 IPLUS)
ENDIF
L0MSKZ=L0MSKZ0 ! restore
ENDIF
C... CAREA is surface area in affected cell.
IF(CAREA>0.0) THEN ! found surface
IF(INDVAR==P1.OR.INDVAR==P2.OR.INDVAR==R1.OR.
1 INDVAR==R2) THEN
C... mass source is density*surface_area_in_cell*velocity_normal_to_surface
IF(QEQ(TOTVEL,-999.)) THEN ! velocity vector has been set
VELMASS=DOT(VELIN,NORML) ! resolve velocity normal to surf
ELSE ! velocity normal to surface has been set
VELMASS=TOTVEL ! velocity is always normal to surface
ENDIF
C... for +ve density, use in-cell density for -ve velocity
C -ve density signals fixed mass flow so take abs
C density=GRND signals fixed volume outflow use in-cell
IF(RHOIN>0.0) THEN ! fixed flow
IF(VELMASS>=0.0) THEN ! inflow
F(L0VAL+I)=RHOIN*ARAT*CAREA*VELMASS
ELSE ! outflow, use in-cell density
F(L0VAL+I)=F(L0DEN+I)*ARAT*CAREA*VELMASS
ENDIF
ELSEIF(QEQ(RHOIN,GRND)) THEN ! fixed volume out flow
F(L0VAL+I)=F(L0DEN+I)*ARAT*CAREA*VELMASS ! - use in-cell
ELSE ! fixed mass flow - use inlet density
F(L0VAL+I)=ABS(RHOIN)*ARAT*CAREA*VELMASS
ENDIF
if(l0vall>0) f(l0vall+i)=f(l0val+i)
ELSEIF(INDVAR==ISURN.OR.INDVAR==LBSRN2) THEN
C... volume source for SURN is surface_area_in_cell*velocity_normal_to_surface
IF(QEQ(TOTVEL,-999.)) THEN ! velocity vector has been set
VELMASS=DOT(VELIN,NORML) ! resolve velocity normal to surf
ELSE ! velocity normal to surface has been set
VELMASS=TOTVEL ! velocity is always normal to surface
ENDIF
F(L0VAL+I)=ARAT*CAREA*VELMASS*SURNEXT
ELSE ! other sources
IF(INDVAR>=U1.AND.INDVAR<=W2) THEN ! momentum sources
C... VAL set to velocity component in relevant direction
IF(CARTES.OR.INDVAR>V2) THEN ! Cartesian or Polar W1/W2
IF(QEQ(TOTVEL,-999.)) THEN ! velocity vector was set
F(L0VAL+I)=VEL
ELSE ! velocity normal to surface has been set
F(L0VAL+I)=TOTVEL*NORML(IDIR) ! resolve into grid dir
ENDIF
ELSE ! Polar U1-V2
IF(QEQ(TOTVEL,-999.)) THEN ! velocity vector was set
UC=VELIN(1) ! set cartesian components
VC=VELIN(2)
ELSE
UC=NORML(1)*TOTVEL ! cart comp is normal*velocity
VC=NORML(2)*TOTVEL
ENDIF
IF(INDVAR0) f(l0vall+i)=f(l0val+i)
ELSEIF(INDVAR==KE) THEN ! KE
F(L0VAL+I)=F(L0KEIN(IDMN)+IPT)
ELSEIF(INDVAR==EP.OR.INDVAR==LBET) THEN ! EP, ET
F(L0VAL+I)=F(L0EPIN(IDMN)+IPT)
ELSEIF(INDVAR==LBOMEG) THEN ! OMEG
F(L0VAL+I)=F(L0EPIN(IDMN)+IPT)/(CMUCD*
1 F(L0KEIN(IDMN)+IPT))
ELSEIF(INDVAR==LBKP) THEN ! KP
F(L0VAL+I)=F(L0KEIN(IDMN)+IPT)/(1+.25)
ELSEIF(INDVAR==LBKT) THEN ! KT
GKTIN=0.25*F(L0KEIN(IDMN)+IPT)/(1+.25)
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDIF
ELSEIF(IGR==13.AND.ISC==15) THEN
C------------------- SECTION 15 ------------------- value = GRND3
C
C Fixed flow rate at Polar INLET
C -------------------------------
C
C The patch can be anything. The COVALs should be:
C COVAL(GXIn, P1, FIXFLU, GRND3)
C COVAL(GXIn, vel, ONLYMS, GRND3)
C
C Inputs:
C SPEDAT(SET,object_name,INLET,C,GXI) - identifies this as an Angled-in patch
C SPEDAT(SET,object_name,UIN,R,U_component)
C SPEDAT(SET,object_name,VIN,R,V_component)
C SPEDAT(SET,object_name,WIN,R,W_component)
C These set the cartesian components of the inlet vector.
C SPEDAT(SET,object_name,VOLFLOW,R,volumetric_flow_rate) - this sets
C the volumetric flow rate. The velocity is deduced from flowrate
C divided by total inlet area.
C SPEDAT(SET,object_name,DENIN,R,inlet_density) - this sets the inlet
C density. if not set, RHO1 is used. If that is GRNDn, 1.0 is used.
C For fixed outflows, the inlet density is stored as -1, to flag the
C use of the in-cell density.
C
IPT=NINT(F(L0GXI(IDMN)+IPAT))
IF(IPT>0) THEN ! current patch is for ANGLED-IN
L0PAT=L0PATNO(IDMN)+(IZ-1)*NXNY ! index for patch-number store
L0VAL=L0F(VAL) ! index for VALue
IF(FLUID1(INDVAR)) THEN
RHOIN=F(L0DENIN(IDMN) +IPT) ! get density at inlet
L0VIN=L0VELI(IDMN)
L0DEN=L0F(DEN1)
ELSE
RHOIN=F(L0DENIN2(IDMN) +IPT) ! get density at inlet
L0VIN=L0VELI2(IDMN)
L0DEN=L0F(DEN2)
ENDIF
CALL GETPAT(IPAT,IREG,TYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
ARAT=F(L0ARATI(IDMN)+IPT) ! area ratio
VELIN(1)=F(L0VIN+(IPT-1)*3+1) ! get U component at inlet
VELIN(2)=F(L0VIN+(IPT-1)*3+2) ! get V component at inlet
VELIN(3)=F(L0VIN+(IPT-1)*3+3) ! get W component at inlet
IF(INDVAR==U1.OR.INDVAR==U2) THEN ! set direction vector
CALL VECTOR(DIR,1.,0.,0.)
ELSEIF(INDVAR==V1.OR.INDVAR==V2) THEN
CALL VECTOR(DIR,0.,1.,0.)
ELSEIF(INDVAR==W1.OR.INDVAR==W2) THEN
CALL VECTOR(DIR,0.,0.,1.)
ENDIF
VEL=DOT(VELIN,DIR) ! get velocity magnitude in velocity direction
IPW=0
L0XU=L0F(XU2D)
L0XG=L0F(XG2D)
IF(IURVAL/=0) L0RAD=L0F(RG2D)
DO IX=IXF,IXL ! loop over all cells in bounding box
DO IY=IYF,IYL
I=(IX-1)*NY+IY
IPW=IPW+1
F(L0VAL+I)=0.0 ! initialise VAL to 0.0
IF(.NOT.SLD(I))THEN
IJKDM=(IX-1)*NY+IY+(IZ-1)*NXNY ! cell index in 3D
IF(INDVAR==P1.OR.INDVAR==P2.OR.INDVAR==R1.OR.
1 INDVAR==R2.OR.INDVAR==ISURN.OR.INDVAR==LBSRN2) THEN
C... mass source is density*velocity_normal_to_surface
IF(QEQ(TYP,2.)) THEN ! East
VELMASS=-VELIN(1)*COS(F(L0XG+I))+
1 VELIN(2)*SIN(F(L0XG+I))
ELSEIF(QEQ(TYP,3.)) THEN ! West
VELMASS= VELIN(1)*COS(F(L0XG+I))-
1 VELIN(2)*SIN(F(L0XG+I))
ELSEIF(QEQ(TYP,4.)) THEN ! North
VELMASS=-VELIN(1)*SIN(F(L0XG+I))-
1 VELIN(2)*COS(F(L0XG+I))
ELSEIF(QEQ(TYP,5.)) THEN ! South
VELMASS= VELIN(1)*SIN(F(L0XG+I))+
1 VELIN(2)*COS(F(L0XG+I))
ELSEIF(QEQ(TYP,6.)) THEN ! High
VELMASS=-VELIN(3)
ELSE ! Low
VELMASS= VELIN(3)
ENDIF
!!! IF(INDVAR/=ISURN) THEN
!!! IF(VELMASS>0.0) THEN ! fixed in flow
!!! F(L0VAL+I)=RHOIN*ARAT*VELMASS
!!! ELSE ! fixed out flow - use in-cell
!!! F(L0VAL+I)=F(L0DEN+I)*ARAT*VELMASS
!!! ENDIF
!!! ELSE
!!!C... volume source for SURN is velocity_normal_to_surface
!!! F(L0VAL+I)=ARAT*VELMASS
!!! ENDIF
IF(INDVAR==ISURN.OR.INDVAR==LBSRN2) THEN
C... volume source for SURN is velocity_normal_to_surface
F(L0VAL+I)=ARAT*VELMASS
ELSE
IF(VELMASS>0.0) THEN ! fixed in flow
F(L0VAL+I)=RHOIN*ARAT*VELMASS
ELSE ! fixed out flow - use in-cell
F(L0VAL+I)=F(L0DEN+I)*ARAT*VELMASS
ENDIF
ENDIF
ELSE ! other sources
IF(INDVAR>=U1.AND.INDVAR<=W2) THEN ! momentum sources
C... VAL set to velocity component in relevant direction
IF(INDVAR>V2) THEN ! Polar W1/W2
F(L0VAL+I)=VEL
ELSE ! Polar U1-V2
IF(INDVAR0) THEN ! link to next
ISOR=IR; ITAR=II
ELSE ! link to previous
ISOR=IABS(IR); ITAR=II
ENDIF
IF(LPAR) THEN ! switch ISOR from global to local
PATNM=GD_NAMPAT(ISOR) ! global patch name
IF(PATNM(1:1)=='^') PATNM=GD_NAMPAT(ISOR)(2:)
ISOR = INDPATCH(PATNM) ! local patch index from global name
ENDIF
IF(ISOR<0) THEN ! No patch in current SubDomain
SMASS1 = 0.0
ELSE ! Patch used, look CoVal(4)
CALL GETSO(ISOR,9,SMASS1) ! get phase 1 mass source
ENDIF
IF(LPAR) CALL GLSUM(SMASS1) ! sum over all nodes
RHOIN=RHO1
IF(RHO1<0.0.AND.QNE(RHO1,GRND5).AND.ITEM1==0) THEN
CALL SET_ERR(573,
1 'Linked ANGLED-IN cannot use current density formula',1)
ENDIF
IF(.NOT.ONEPHS) THEN
IF(ISOR<0) THEN ! No patch in current SubDomain
SMASS2 = 0.0
ELSE ! Patch used, look CoVal(4)
CALL GETSO(ISOR,10,SMASS2)
ENDIF
IF(LPAR) CALL GLSUM(SMASS2) ! sum over all nodes
RHOIN2=RHO2
IF(RHO2<0.0.AND.QNE(RHO2,GRND5).AND.ITEM2==0) THEN
CALL SET_ERR(573,
1 'Linked ANGLED-IN cannot use current density formula',1)
ENDIF
ENDIF
IF(ITAR>0) ITR=NINT(F(L0GXI(IDMN)+ITAR)) ! sequence number of target
DO II=1,NSOL
IPHI=MSL(II); IF(IPHI<14) CYCLE
IF(ISOR<0)THEN ! source patch not on this processor
GSOR=0.0; GCO=0.0
ELSE
CALL GETCV(ISOR,IPHI,GCO,GVAL) ! get C V from COVAL
IF(QEQ(GCO,-999.)) THEN ! no source was set so skip
IF(LPAR) THEN
GSOR=0.0 ! set zero value for parallel sum
ELSE
CYCLE ! skip everything for sequential
ENDIF
ELSE
CALL GETSO(ISOR,IPHI,GSOR) ! get source of variable in 'donor'
ENDIF
ENDIF
IF(LPAR) THEN
CALL GLSUM2(GSOR,GCO) ! sum source and CO over all processors
IF(GCO<-998.0) CYCLE ! if CO is < -998, must be -999 somewhere so skip
ENDIF
IF(ITAR<0) CYCLE ! target not on this processor
IF(ONEPHS.OR.FLUID1(IPHI)) THEN
SMASS=SMASS1
ELSE
SMASS=SMASS2
ENDIF
IF(NAME(IPHI)=='TEM1') THEN ! Phase 1 temperature
IF(CP1<0.0) CALL SET_ERR(572,
1 'Linked ANGLED-IN cannot use variable specific heat',1)
IQTAD1=F(L0IADQ(IDMN)+ITR) ! source-type flag
IF(IQTAD1==0) THEN ! additional heat
QADD=F(L0ADDQ(IDMN)+ITR)
GAVE=(GSOR-QADD)/SMASS/CP1-TEMP0 ! subtract as SMASS is -ve!
ELSEIF(IQTAD1==1) THEN ! set exit temperature
GAVE=F(L0ADDQ(IDMN)+ITR)
ELSEIF(IQTAD1==2) THEN ! additional temperature
GAVE=GSOR/SMASS/CP1-TEMP0+F(L0ADDQ(IDMN)+ITR)
ENDIF
TMIN=F(L0ADMIN(IDMN,IPHI)+ITR)
TMAX=F(L0ADMAX(IDMN,IPHI)+ITR)
GAVE=MAX(TMIN,MIN(GAVE,TMAX))
T1AVE=GAVE
ELSEIF(NAME(IPHI)=='TEM2') THEN ! Phase 2 temperaure
IF(CP2<0.0) CALL SET_ERR(572,
1 'Linked ANGLED-IN cannot use variable specific heat',1)
IQTAD1=F(L0IADQ(IDMN)+ITR)
IF(IQTAD1==0) THEN ! additional heat
QADD=F(L0ADDQ(IDMN)+ITR)
GAVE=(GSOR-QADD)/SMASS/CP2-TEMP0
ELSEIF(IQTAD1==1) THEN ! set exit temperature
GAVE=F(L0ADDQ(IDMN)+ITR)
ELSEIF(IQTAD1==2) THEN ! additional temperature
GAVE=GSOR/SMASS/CP2-TEMP0+F(L0ADDQ(IDMN)+ITR)
ENDIF
TMIN=F(L0ADMIN(IDMN,IPHI)+ITR)
TMAX=F(L0ADMAX(IDMN,IPHI)+ITR)
GAVE=MAX(TMIN,MIN(GAVE,TMAX))
T2AVE=GAVE
ELSEIF(IPHI==14.OR.IPHI==15) THEN ! enthalpy
QADD=F(L0ADDQ(IDMN)+ITR)
GAVE=(GSOR-QADD)/SMASS
TMIN=F(L0ADMIN(IDMN,IPHI)+ITR)
TMAX=F(L0ADMAX(IDMN,IPHI)+ITR)
GAVE=MAX(TMIN,MIN(GAVE,TMAX))
ELSE ! other scalars
ICAD=F(L0IADC(IDMN,IPHI)+ITR) ! source-type flag
IF(NAME(IPHI)=='MH2O') THEN ! deal with Humidity equation MH2O
IHUNIT=F(L0HUNIT(IDMN)+ITR) ! get flag for humidity source units
CADD=F(L0ADDC(IDMN,IPHI)+ITR) ! get humidity source value
IF(IHUNIT>0.AND.ICAD<3) THEN ! not mass fraction andnot % reduction
HUMIN=CADD
IF(IHUNIT==1) THEN ! convert from humidity ratio
HUMIN=1.E-3*HUMIN
ELSEIF(IHUNIT==2) THEN ! convert from relative humidity
CALL GETSO(ITAR,ITEM1,GTSOR) ! energy source at this patch
CALL GETSO(ITAR,P1,GMSOR) ! mass source at this patch
TIN=GTSOR/(GMSOR+TINY)/CP1-TEMP0 ! outlet temperature (may be 1 sweep behind)
PVAP = 1.E-2*HUMIN*PVH2O(TIN)
GWRAT = 18.015/28.96
HUMIN = GWRAT*PVAP/(PRESS0-PVAP)
ENDIF
CADD=HUMIN/(1.+HUMIN) ! convert to mass fraction
ENDIF
IF(ICAD==0) THEN ! extra humidity source (flux)
GAVE=(GSOR-CADD)/SMASS ! subtract as SMASS is -ve!
ELSEIF(ICAD==1) THEN ! Set exit humidity
GAVE=CADD
ELSEIF(ICAD==2) THEN ! Set additional exit humidity
GAVE=GSOR/SMASS+CADD
ELSEIF(ICAD==3) THEN ! Set % humidity reduction - treat as % mass fraction reducti
RFAC=1.-(F(L0ADDC(IDMN,IPHI)+ITR)/100.)
GAVE=(GSOR/SMASS)*RFAC ! Reduce by %
ENDIF
CMIN=F(L0ADMIN(IDMN,IPHI)+ITR)
CMAX=F(L0ADMAX(IDMN,IPHI)+ITR)
GAVE=MAX(CMIN,MIN(GAVE,CMAX))
ELSE ! All other scalars not handled above
IF(ICAD==0) THEN ! additional c source
CADD=F(L0ADDC(IDMN,IPHI)+ITR)
GAVE=(GSOR-CADD)/SMASS ! subtract as SMASS is -ve!
ELSEIF(ICAD==1) THEN ! set exit C
GAVE=F(L0ADDC(IDMN,IPHI)+ITR)
ELSEIF(ICAD==2) THEN ! additional C value
GAVE=GSOR/SMASS+F(L0ADDC(IDMN,IPHI)+ITR)
ELSEIF(ICAD==3) THEN ! set % reduction
RFAC=1.-(F(L0ADDC(IDMN,IPHI)+ITR)/100.)
GAVE=(GSOR/SMASS)*RFAC ! Reduce by %
ENDIF
CMIN=F(L0ADMIN(IDMN,IPHI)+ITR)
CMAX=F(L0ADMAX(IDMN,IPHI)+ITR)
GAVE=MAX(CMIN,MIN(GAVE,CMAX))
IF(IPHI==C1) THEN
C1AVE=GAVE
ELSEIF(IPHI==C2) THEN
C2AVE=GAVE
ENDIF
ENDIF
ENDIF
CALL CORCOVAL(ITAR,IPHI,GCO,GAVE) ! set value back
ENDDO
IF(ITAR<0) CYCLE ! target not on this processor
IF(QEQ(RHO1,GRND5)) THEN
IF(EQZ(RHO1A)) THEN
RHOIN=RHO1B*PRESS0/(T1AVE+TEMP0) ! inlet density
ELSE
SPEVOL=RHO1A+(RHO1B-RHO1A)*C1AVE
RHOIN=PRESS0/(SPEVOL*(T1AVE+TEMP0))
ENDIF
ENDIF
IF(QEQ(RHO2,GRND5)) THEN
IF(EQZ(RHO2A)) THEN
RHOIN2=RHO2B*PRESS0/(T2AVE+TEMP0)
ELSE
SPEVOL=RHO2A+(RHO2B-RHO2A)*C2AVE
RHOIN2=PRESS0/(SPEVOL*(T2AVE+TEMP0))
ENDIF
ENDIF
ARAT=F(L0ARATI(IDMN)+ITR) ! area ratio
TOTAREA=F(L0TOTAREA(IDMN)+ITR) ! area
C... Deduce discharge velocity for target as mass/(area*density)
F(L0TOTVEL(IDMN)+ITR)=ABS(SMASS1)/(ARAT*RHOIN*TOTAREA+TINY)
C... for mass flow rate, set density to -ve, to flag mass flow
F(L0DENIN(IDMN)+ITR)=-RHOIN
IF(.NOT.ONEPHS) THEN
F(L0TOTVEL2(IDMN)+ITR)=ABS(SMASS2)/(ARAT*RHOIN2*TOTAREA+
1 TINY)
F(L0DENIN2(IDMN)+ITR)=-RHOIN2
ENDIF
IF(GRN(ENUT)) THEN
C... Inlet values for Turbulence
TURBIN=F(L0TURBIN(IDMN)+ITR)
IF(QNE(TURBIN,-999.0)) THEN ! intensity was set, recompute KE & EP inlet
IF(SOLVE(KE)) THEN
GKEIN=(TURBIN*F(L0TOTVEL(IDMN)+ITR))**2
F(L0KEIN(IDMN)+ITR)=GKEIN
ENDIF
IF(SOLVE(EP).OR.SOLVE(LBOMEG).OR.SOLVE(LBET).OR.
1 SOLVE(LBKP) .OR.SOLVE(LBKT)) THEN
RHYD=(F(L0TOTAREA(IDMN)+ITR)/(4.*ATAN(1.)))**0.5
F(L0EPIN(IDMN)+ITR)=CMUCD**0.75*GKEIN**1.5/(0.1*RHYD)
ENDIF
ENDIF
ENDIF ! end enut==grnd
ENDIF ! end 'if linked'
ENDDO ! end loop oer angled-ins
ENDIF
C****************************************************************
END
C---------------------------------------------------------------------
SUBROUTINE GET_SURFACE(IX,IY,IZ,IBLK,L0PAT,CAREA,NORML,IPLUS)
C... This subroutine returns the surface area and surface normal for
C a cell lying on the intersection of a blockage and the current
C patch.
C Inputs:
C IX,IY,IZ - cell index
C IBLK - patch number of blockage. If 0, any blockage will be used
C L0PAT - address of patch-number store
C outputs:
C CAREA - surface area
C NORML(3) - surface normal unit vector
C IPLUS - 1 if surface area is at IZ+1 for W1
C
COMMON/ISG/ISGF1(61),ISG62,ISGF2(38) ! isg62==1 for sparsol
REAL NORML(3)
C
IF(ISG62==0) THEN
CALL GET_SURFACE_P(IX,IY,IZ,IBLK,L0PAT,CAREA,NORML,IPLUS)
ELSE
CALL GET_SURFACE_S(IX,IY,IZ,IBLK,L0PAT,CAREA,NORML,IPLUS)
ENDIF
END
C---------------------------------------------------------------------
SUBROUTINE GET_SURFACE_P(IX,IY,IZ,IBLK,L0PAT,CAREA,NORML,IPLUS)
C... This subroutine returns the surface area and surface normal for
C a cell lying on the intersection of a blockage and the current
C patch. For Parsol and Xparsol.
C Inputs:
C IX,IY,IZ - cell index
C IBLK - patch number of blockage. If 0, any blockage will be used
C L0PAT - address of patch-number store
C outputs:
C CAREA - surface area
C NORML(3) - surface normal unit vector
C IPLUS - 1 if surface area is at IZ+1 for W1
C
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
COMMON/LBFFV/P1,P2,U1,U2,V1,V2,W1,W2,R1,R2,RS,KE,EP,H1,H2,
1C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,
1C18,C19,C20,C21,C22,C23,C24,C25,C26,C27,C28,C29,C30,C31,C32,
1C33,C34,C35
INTEGER P1,P2,U1,U2,V1,V2,W1,W2,R1,R2,RS,KE,EP,H1,H2,C1,C2,C3,
1C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,C18,C19,C20,
1C21,C22,C23,C24,C25,C26,C27,C28,C29,C30,C31,C32,C33,C34,C35
INCLUDE 'd_earth/pbcstr'
COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
COMMON /INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT,
1 NDOMMX
COMMON /GEODMN0/ I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/F0/KFIL1(17),KZZG, KFIL2(17), KXYXG,KXYXU,KXYYG, KFIL3(266)
COMMON/ISG/ISGF1(61),ISG62,ISGF2(38) ! isg62==1 for sparsol
REAL NORML(3),norm0(3),norm1(3)
LOGICAL DOIT,EQZ,SLD
IPLUS=0
I=(IX-1)*NY+IY ! cell index in plane
IJKDM=(IX-1)*NY+IY+(IZ-1)*NXNY ! cell index in 3D
CAREA=0.0 ! initialise surface area
DOIT=.FALSE. ! initialise 'proceed' flag
IPBC=0
IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
IF(IPBC>0) THEN ! dealing with cut cell
IF(IBLK==0) THEN ! no blockage specified - any will do
DOIT=F(L0PAT+I)>0.0
ELSE
DOIT=NINT(F(L0PAT+I))==IBLK
ENDIF
IF(DOIT) THEN ! this is the right cell!
ICUTIDX = ISHPB(IDMN,IPBC)
CAREA=F(PBC_CTAR+ICUTIDX) ! surface area of cut cell
CALL VECTOR(NORML,F(PBC_CTVEC(1)+ICUTIDX),
1 F(PBC_CTVEC(2)+ICUTIDX),F(PBC_CTVEC(3)+ICUTIDX))
CALL NORM(NORML) ! unit-normal to cut surface
ENDIF
ELSE ! dealing with un-cut cell
INE=ITWO(I+NY, I,IX 1) ! - West
INN=ITWO(I+1, I,IY 1) ! - South
INH=ITWO(I+NXNY,I,IZ 1) ! - Low
IF(F(L0PAT+INE)>0.0) THEN ! solid to east
IF(IBLK==0) THEN ! any blockage will do
DOIT=.TRUE.
ELSE ! only next to blockage IBLK
DOIT=NINT(F(L0PAT+INE))==IBLK
ENDIF
IF(DOIT) THEN
IJKDME=IJKDM+NY ! neighbour index
IPBCE = IPBSEQ(IDMN,IJKDME) ! neighbour cut-cell index
IF(IPBCE>0) THEN ! neighbour is cut-cell
ICUTIDXE = ISHPB(IDMN,IPBCE)
AREA=F(PBC_DAKAR(4)+ICUTIDXE) ! surface area for cut-cell
ELSE ! neighbour is fully blocked
AREA=F(I3DAEX+IJKDM) ! surface area from cell area
ENDIF
IF(AREA>CAREA) THEN
CAREA=AREA
IF(CARTES) THEN
CALL VECTOR(NORML,-1.,0.,0.) ! set unit vector normal to surface
ELSE
XC=-COS(F(KXYXU+I))
YC= SIN(F(KXYXU+I))
CALL VECTOR(NORML,XC,YC,0.0)
ENDIF
ENDIF
ENDIF
ENDIF
IF(F(L0PAT+INW)>0.0) THEN ! solid to west
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INW))==IBLK
ENDIF
IF(DOIT) THEN
IJKDMW=IJKDM-NY
IPBCW = IPBSEQ(IDMN,IJKDMW)
IF(IPBCW>0) THEN
ICUTIDXW = ISHPB(IDMN,IPBCW)
AREA=F(PBC_DAKAR(3)+ICUTIDXW)
ELSE
AREA=F(I3DAEX+IJKDMW)
ENDIF
IF(AREA>CAREA) THEN
CAREA=AREA
IF(CARTES) THEN
CALL VECTOR(NORML,1.,0.,0.)
ELSE
XC= COS(F(KXYXU+I))
YC= -SIN(F(KXYXU+I))
CALL VECTOR(NORML,XC,YC,0.0)
ENDIF
ENDIF
ENDIF
ENDIF
IF(F(L0PAT+INN)>0.0) THEN ! solid to north
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INN))==IBLK
ENDIF
IF(DOIT) THEN
IJKDMN=IJKDM+1
IPBCN = IPBSEQ(IDMN,IJKDMN)
IF(IPBCN>0) THEN
ICUTIDXN = ISHPB(IDMN,IPBCN)
AREA=F(PBC_DAKAR(2)+ICUTIDXN)
ELSE
AREA=F(I3DANY+IJKDM)
ENDIF
IF(AREA>CAREA) THEN
CAREA=AREA
IF(CARTES) THEN
CALL VECTOR(NORML,0.,-1.,0.)
ELSE
XC=-SIN(F(KXYXG+I))
YC=-COS(F(KXYXG+I))
CALL VECTOR(NORML,XC,YC,0.0)
ENDIF
ENDIF
ENDIF
ENDIF
IF(F(L0PAT+INS)>0.0) THEN ! solid to south
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INS))==IBLK
ENDIF
IF(DOIT) THEN
IJKDMS=IJKDM-1
IPBCS = IPBSEQ(IDMN,IJKDMS)
IF(IPBCS>0) THEN
ICUTIDXS = ISHPB(IDMN,IPBCS)
AREA=F(PBC_DAKAR(1)+ICUTIDXS)
ELSE
AREA=F(I3DANY+IJKDMS)
ENDIF
IF(AREA>CAREA) THEN
CAREA=AREA
IF(CARTES) THEN
CALL VECTOR(NORML,0.,1.,0.)
ELSE
XC=SIN(F(KXYXG+I-1))
YC=COS(F(KXYXG+I-1))
CALL VECTOR(NORML,XC,YC,0.0)
ENDIF
ENDIF
ENDIF
ENDIF
IF(F(L0PAT+INH)>0.0) THEN ! solid to high
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INH))==IBLK
ENDIF
IF(DOIT) THEN
IJKDMH=IJKDM+NXNY
IPBCH = IPBSEQ(IDMN,IJKDMH)
IF(IPBCH>0) THEN
ICUTIDXH = ISHPB(IDMN,IPBCH)
AREA=F(PBC_DAKAR(6)+ICUTIDXH)
ELSE
AREA=F(I3DAHZ+IJKDM)
ENDIF
IF(AREA>CAREA) THEN
CAREA=AREA
CALL VECTOR(NORML,0.,0.,-1.)
ENDIF
ENDIF
ENDIF
IF(F(L0PAT+INL)>0.0) THEN ! solid to low
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INL))==IBLK
ENDIF
IF(DOIT) THEN
IJKDML=IJKDM-NXNY
IPBCL = IPBSEQ(IDMN,IJKDML)
IF(IPBCL>0) THEN
ICUTIDXL = ISHPB(IDMN,IPBCL)
AREA=F(PBC_DAKAR(5)+ICUTIDXL)
ELSE
AREA=F(I3DAHZ+IJKDML)
ENDIF
IF(AREA>CAREA) THEN
CAREA=AREA
CALL VECTOR(NORML,0.,0.,1.)
ENDIF
ENDIF
ENDIF
ENDIF
IF((INDVAR==W1.OR.INDVAR==W2).AND.EQZ(CAREA).AND.IZ0) THEN ! dealing with cut cell
IF(IBLK==0) THEN
DOIT=F(L0PAT1+I)>0.0
ELSE
DOIT=NINT(F(L0PAT1+I))==IBLK
ENDIF
IF(DOIT) THEN
ICUTIDX = ISHPB(IDMN,IPBC)
CAREA=F(PBC_CTAR+ICUTIDX) ! surface area of cut cell
CALL VECTOR(NORML,F(PBC_CTVEC(1)+ICUTIDX),
1 F(PBC_CTVEC(2)+ICUTIDX),F(PBC_CTVEC(3)+ICUTIDX))
CALL NORM(NORML) ! unit-normal to cut surface
ENDIF
ELSEIF(IZ0.0) THEN ! solid to high
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT1+INH))==IBLK
ENDIF
IF(DOIT) THEN
CALL VECTOR(NORML,0.,0.,-1.)
IJKDMH=IJKDM+NXNY
IPBCH = IPBSEQ(IDMN,IJKDMH)
IF(IPBCH>0) THEN
ICUTIDXH = ISHPB(IDMN,IPBCH)
CAREA=F(PBC_DAKAR(6)+ICUTIDXH)
ELSE
CAREA=F(I3DAHZ+IJKDMH)
ENDIF
ENDIF
ENDIF
ENDIF
IF(CAREA>0) IPLUS=1
ENDIF
END
C---------------------------------------------------------------------
SUBROUTINE GET_SURFACE_S(IX,IY,IZ,IBLK,L0PAT,CAREA,NORML,IPLUS)
C... This subroutine returns the surface area and surface normal for
C a cell lying on the intersection of a blockage and the current
C patch. For Sparsol.
C Inputs:
C IX,IY,IZ - cell index
C IBLK - sequence number of blockage. If 0, any blockage will be used
C L0PAT - address of patch-number store
C outputs:
C CAREA - surface area
C NORML(3) - surface normal unit vector
C IPLUS - 1 if surface area is at IZ+1 for W1
C
use cut_cell_data
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
COMMON/LBFFV/P1,P2,U1,U2,V1,V2,W1,W2,R1,R2,RS,KE,EP,H1,H2,
1C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,
1C18,C19,C20,C21,C22,C23,C24,C25,C26,C27,C28,C29,C30,C31,C32,
1C33,C34,C35
INTEGER P1,P2,U1,U2,V1,V2,W1,W2,R1,R2,RS,KE,EP,H1,H2,C1,C2,C3,
1C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,C18,C19,C20,
1C21,C22,C23,C24,C25,C26,C27,C28,C29,C30,C31,C32,C33,C34,C35
INCLUDE 'd_earth/pbcstr'
COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
COMMON /INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT,
1 NDOMMX
COMMON /GEODMN0/ I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
! COMMON/F0/KFIL1(17),KZZG, KFIL2(17), KXYXG,KXYXU,KXYYG, KFIL3(266)
COMMON/F0/KXDX,KXXG,KXXU,KXDXG,KXS,KYDY,KYYG,KYR,KYR2,KYRV,KYRV2,
1K12CNN,K12DFN,KYYV,KYDYG,KYS,KZDZ,KZZG,KZP,KZPP,KZAP,KZAH,KZSU,
1KZAL,KZZW,KZDZG,KZWG,KZBLK,KZDZO,KXYDX,KXYDY,KTZPRF,KXYR,KD11,
1KXYRV,KXYXG,KXYXU,KXYYG,KXYYV,KXYDXG,KXYDYG,KZALF,K3PHI,KZAPF,
1KZAHF,KZSUF,KABS,KOR1,K3BB,KCM1,KOR2,KMULM,KMIN1,KMOUT1,KM1,
1K12CON,KGAM,KGAMH,KCLA,KMU,KMUH,K12SOR,KCLA2,KST1,KAE,KAW,KAN,KAS,
1KAH,KAL,KAP,KSU,KARL,IPBIGD,KCM2,KMIN2,KMOUT2,KM2,KLW1AD,KLW2AD,
1KGEN1,KGEN2,K12CNH,IPAWD,IPASD,IPALD,KPLOT,KST2,KAP2,KSU2,KMDOT,
1KS12,IPAED,IPAND,KD12,KMULMH,KD1,KD2,KD3,KD4,KD5,KD6,KD7,KD8,KD9,
1KD10,KDSTAG,KCP1L,KR1H,KZXCY,KCP2L,K12DFH,KU,KU2,KV,KV2,KW,KW2,
1KDU,KDU2,KDV,KDV2,KDW,KDW2,KCE,KCW,KCN,KCS,KCH,KCL,KCE2,KCW2,KCN2,
1KCS2,KCH2,KCL2,KCP1,KRHO1,KRHO1H,KCP2,KRHO2,KRHO2H,IPAHD,KR1,IVOL,
1IAEX,IANY,IAHZ,KXYGR(10),KXYAE,KXYAW,KXYAN,KXYAS,KXYAH,KXYAL,
1KXYVOL,KXYDZ,KXYDZG,KXYDZL,KCP1H,KT1,KT1H,KCP2H,KT2,KT2H,K12CNE,
1KL1,KL1H,K12DFE,KL2,KL2H,IPAE,IPAN,IPAH,IPAP,IPSU,IQDU,IQDU2,IQDV,
1IQDV2,IQDW,IQDW2,IPAS,IPAW,IPAL,KCEA,KCNA,KCHA,KCEA2,KCNA2,KCHA2,
1KX1D,KX2D,KXD,KXAPF,KXAEF,KXAWF,KXSUF,KXPP,KY1D,KY2D,KYD,KYAPF,
1KYANF,KYASF,KYSUF,KYPP,KZ1D,KZ2D,KZD,KYZAW,KYZAE,KXYEA(20),KACO,
1KUCO,KLCO,KCR,KRS,KZYCYP,KDRE,KDRN,KDRH,KRHOM,KRHOMH,KCP1O,KCP2O,
1KGMLM,KGMLMH,IVOLH,KATMP,I1E,I2E,I1N,I2N,I1H,I2H,IFACE,IF1(39)
COMMON/STOGEO/STORGEO
LOGICAL STORGEO
COMMON/ISG/ISGF1(61),ISG62,ISGF2(38) ! isg62==1 for sparsol
COMMON/L0VRPAT/L0IFLP ! for the IFLP index
COMMON/L0PVCM/IDUM1(8),L0FLP2,NPTINF
REAL NORML(3),norm0(3),norm1(3)
INTEGER MAXDIR(1)
LOGICAL DOIT,EQZ,SLD,DOIT2
IPLUS=0
I=(IX-1)*NY+IY ! cell index in plane
IJKDM=(IX-1)*NY+IY+(IZ-1)*NXNY ! cell index in 3D
CAREA=0.0; NORML=0.0 ! initialise surface area and normal
DOIT=.FALSE. ! initialise 'proceed' flag
DOIT2=ISWEEP==FSWEEP.AND.ISTEP==FSTEP.AND.(INDVARW2)
ICUT=0
IF(ASSOCIATED(CUTCELL)) ICUT = CUTCELL%ICUTM(IJKDM)
IF(ICUT>0) THEN ! dealing with cut cell
IF(IBLK==0) THEN ! no blockage specified - any will do
DOIT=F(L0PAT+I)>0.0.OR.ICUT>0
ELSE
IFLP=F(L0IFLP+IPAT) ! number of facet object connected with ipat
IBLKD=NINT(ABS(F(L0FLP2+(IFLP-1)*6+1)))
DOIT=CUTCELL%IOBJ(ICUT)==IBLKD
ENDIF
IF(DOIT) THEN ! this is the right cell!
CAREA=CUTCELL%AREA(ICUT) ! area of cut surface
CALL VECTOR(NORML,CUTCELL%NORML(1,ICUT),
1 CUTCELL%NORML(2,ICUT),CUTCELL%NORML(3,ICUT))
IF(DOIT2) THEN ! undo Sparsol area adjustment,restore nominal areas.
IF(NX>1) THEN
AREAE=F(KYDY+IY)*F(KZDZ+IZ)
F(I3DAEX+IJKDM)=AREAE
IF(IX>1) F(I3DAEX+IJKDM-NY)=AREAE
ENDIF
IF(NY>1) THEN
AREAN=F(KXDX+IX)*F(KZDZ+IZ)
F(I3DANY+IJKDM)=AREAN
IF(IY>1) F(I3DANY+IJKDM-1)=AREAN
ENDIF
IF(NZ>1) THEN
AREAH=F(KXDX+IX)*F(KYDY+IY)
F(I3DAHZ+IJKDM)=AREAH
IF(IZ>1) F(I3DAHZ+IJKDM-NXNY)=AREAH
ENDIF
ENDIF
ENDIF
CALL NORM(NORML) ! unit-normal to cut surface
ELSE ! dealing with un-cut cell
INE=ITWO(I+NY, I,IX 1) ! - West
INN=ITWO(I+1, I,IY 1) ! - South
INH=ITWO(I+NXNY,I,IZ 1) ! - Low
NORML=0.0; NORM0=0.0
IF(IX0.0) THEN ! solid to east
IF(IBLK==0) THEN ! any blockage will do
DOIT=.TRUE.
ELSE ! only next to blockage IBLK
DOIT=NINT(F(L0PAT+INE))==IBLK
ENDIF
C.... don't set area if neighbour is not solid
IF(.NOT.SLD(INE)) DOIT=.FALSE.
IF(DOIT) THEN
IJKDME=IJKDM+NY ! neighbour index
AREA=F(I3DAEX+IJKDM) ! surface area from cell area
IF(ASSOCIATED(CUTCELL)) THEN
IF(CUTCELL%ICUTM(IJKDM+NY)>0) THEN ! neighbour is cut, but is now filled
! source from neighbour will have been moved
! away from surface to an open cell
AREA=AREA*0.001 ! set small source area to prevent 'wind' along wall
IF(DOIT2) THEN
AREAW=F(KYDY+IY)*F(KZDZ+IZ) ! undo Sparsol area adjustment
F(I3DAEX+IJKDM-NY)=AREAW
ENDIF
ENDIF
ENDIF
CAREA=CAREA+AREA ! add area to total
IF(CARTES) THEN
CALL VECTOR(NORML,-1.,0.,0.) ! set unit vector normal to surface
ELSE
XC=-COS(F(KXYXU+I)); YC= SIN(F(KXYXU+I))
CALL VECTOR(NORML,XC,YC,0.0); CALL NORM(NORML)
ENDIF
CALL ADDVEC(NORM0,NORML,NORM0) ! add normal to total
ENDIF
ENDIF
C
IF(IX>1.AND.F(L0PAT+INW)>0.0) THEN ! solid to west
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INW))==IBLK
ENDIF
C.... don't set area if neighbour is not solid
IF(.NOT.SLD(INW)) DOIT=.FALSE.
IF(DOIT) THEN
IJKDMW=IJKDM-NY
AREA=F(I3DAEX+IJKDMW)
IF(ASSOCIATED(CUTCELL)) THEN
IF(CUTCELL%ICUTM(IJKDM-NY)>0) THEN ! neighbour is cut, but is now filled
AREA=0.001*AREA
IF(DOIT2) THEN
AREAE=F(KYDY+IY)*F(KZDZ+IZ)
F(I3DAEX+IJKDM)=AREAE
ENDIF
ENDIF
ENDIF
CAREA=CAREA+AREA
IF(CARTES) THEN
CALL VECTOR(NORML,1.,0.,0.)
ELSE
XC= COS(F(KXYXU+I)); YC= -SIN(F(KXYXU+I))
CALL VECTOR(NORML,XC,YC,0.0); CALL NORM(NORML)
ENDIF
CALL ADDVEC(NORM0,NORML,NORM0)
ENDIF
ENDIF
C
IF(IY0.0) THEN ! solid to north
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INN))==IBLK
ENDIF
IF(.NOT.SLD(INN)) DOIT=.FALSE.
IF(DOIT) THEN
IJKDMN=IJKDM+1
AREA=F(I3DANY+IJKDM) ! surface area from cell area
IF(ASSOCIATED(CUTCELL)) THEN
IF(CUTCELL%ICUTM(IJKDM+1)>0) THEN
AREA=0.001*AREA
IF(DOIT2) THEN
AREAS=F(KXDX+IX)*F(KZDZ+IZ)
F(I3DANY+IJKDM-1)=AREAS
ENDIF
ENDIF
ENDIF
CAREA=CAREA+AREA
IF(CARTES) THEN
CALL VECTOR(NORML,0.,-1.,0.)
ELSE
XC=-SIN(F(KXYXG+I)); YC=-COS(F(KXYXG+I))
CALL VECTOR(NORML,XC,YC,0.0); CALL NORM(NORML)
ENDIF
CALL ADDVEC(NORM0,NORML,NORM0)
ENDIF
ENDIF
C
IF(IY>1.AND.F(L0PAT+INS)>0.0) THEN ! solid to south
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INS))==IBLK
ENDIF
C.... don't set area if neighbour is not solid
IF(.NOT.SLD(INS)) DOIT=.FALSE.
IF(DOIT) THEN
IJKDMS=IJKDM-1
AREA=F(I3DANY+IJKDMS)
IF(ASSOCIATED(CUTCELL)) THEN
IF(CUTCELL%ICUTM(IJKDM-1)>0) THEN
AREA=0.001*AREA
IF(DOIT2) THEN
AREAN=F(KXDX+IX)*F(KZDZ+IZ)
F(I3DANY+IJKDM)=AREAN
ENDIF
ENDIF
ENDIF
CAREA=CAREA+AREA
IF(CARTES) THEN
CALL VECTOR(NORML,0.,1.,0.)
ELSE
XC=SIN(F(KXYXG+I-1)); YC=COS(F(KXYXG+I-1))
CALL VECTOR(NORML,XC,YC,0.0); CALL NORM(NORML)
ENDIF
CALL ADDVEC(NORM0,NORML,NORM0)
ENDIF
ENDIF
C
IF(IZ0.0) THEN ! solid to high
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INH))==IBLK
ENDIF
C... don't set area if neighbour is not solid
IF(.NOT.SLD(INH)) DOIT=.FALSE.
IF(DOIT) THEN
IJKDMH=IJKDM+NXNY
AREA=F(I3DAHZ+IJKDM)
IF(ASSOCIATED(CUTCELL)) THEN
IF(CUTCELL%ICUTM(IJKDM+NXNY)>0) THEN
AREA=AREA*0.001
IF(DOIT2) THEN
AREAL=F(KXDX+IX)*F(KYDY+IY)
F(I3DAHZ+IJKDM-NXNY)=AREAL
ENDIF
ENDIF
ENDIF
CAREA=CAREA+AREA
CALL VECTOR(NORML,0.,0.,-1.)
CALL ADDVEC(NORM0,NORML,NORM0)
ENDIF
ENDIF
C
IF(IZ>1.AND.F(L0PAT+INL)>0.0) THEN ! solid to low
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT+INL))==IBLK
ENDIF
C.... don't set area if neighbour is not solid
IF(.NOT.SLD(INL)) DOIT=.FALSE.
IF(DOIT) THEN
IJKDML=IJKDM-NXNY
AREA=F(I3DAHZ+IJKDML)
IF(ASSOCIATED(CUTCELL)) THEN
IF(CUTCELL%ICUTM(IJKDM-NXNY)>0) then
AREA=0.001*AREA
IF(DOIT2) THEN
AREAH=F(KXDX+IX)*F(KYDY+IY)
F(I3DAHZ+IJKDM)=AREAH
ENDIF
ENDIF
ENDIF
CAREA=CAREA+AREA
CALL VECTOR(NORML,0.,0.,1.)
CALL ADDVEC(NORM0,NORML,NORM0)
ENDIF
ENDIF
CALL NORM(NORM0); NORML=NORM0
MAXDIR=MAXLOC((/ABS(NORML(1)),ABS(NORML(2)),ABS(NORML(3))/))
CAREA=CAREA*ABS(NORML(MAXDIR(1)))
ENDIF
IF((INDVAR==W1.OR.INDVAR==W2).AND.EQZ(CAREA).AND.IZ0.0) THEN ! solid to high
IF(IBLK==0) THEN
DOIT=.TRUE.
ELSE
DOIT=NINT(F(L0PAT1+INH))==IBLK
ENDIF
IF(DOIT) THEN
CALL VECTOR(NORML,0.,0.,-1.)
IJKDMH=IJKDM+NXNY
CAREA=F(I3DAHZ+IJKDMH)
IF(ASSOCIATED(CUTCELL)) THEN
IF(CUTCELL%ICUTM(IJKDM+NXNY)>0) THEN
CAREA=CAREA*0.001
IF(DOIT2) THEN
AREAL=F(KXDX+IX)*F(KYDY+IY)
F(I3DAHZ+IJKDM-NXNY)=AREAL
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
IF(CAREA>0) IPLUS=1
ENDIF
END
c