c

C.... FILE NAME GXIO.FTN-------------------------------- 240810
      SUBROUTINE GXIO
      INCLUDE '/phoenics/d_includ/patnos'
      INCLUDE '/phoenics/d_includ/patcmn'
      INCLUDE '/phoenics/d_includ/objnam'
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoeclos/d_includ/d_earth/parvar'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      INCLUDE '/phoeclos/d_includ/d_earth/pbcstr'
      INCLUDE '/phoenics/d_includ/parear'
      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 /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
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)
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
      SAVE LPAR
C
C***********************************************************************
      IF(IGR.EQ.1.AND.ISC.EQ.2) THEN
        IF(NUMDMN.GT.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.GT.1) CALL INDXDM
        NGXIN(IDMN)=0
        NGXOUT(IDMN)=0
        NCELL=0
        LPAR=MIMD.AND.NPROC.GT.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.LT.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.EQ.'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.EQ.'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.EQ.'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.EQ.'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).GT.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.GT.0) THEN
            CALL GXMAKE(L0REL(IDMN),NUM,'RELX')  ! relaxation for deduced
            CALL GXMAKE(L0VEL(IDMN),NUM*NCELL,'VELIN') ! initial guess for deduc
          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.GT.0) THEN
              CALL GXMAKE(L0RL2(IDMN),NUM,'RELX2')  ! relaxation for deduced
              CALL GXMAKE(L0VL2(IDMN),NUM*NCELL,'VELIN2') ! initial guess for de
            ENDIF
          ENDIF
        ENDIF
C
        IF(NGXIN(IDMN).GT.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')
          ENDIF
          CALL GXMAKE0(L0TOTAREA(IDMN),NUM,'TOTAREA')
          NUM2=ITWO(ILIM,NUM,LPAR) ! parallel LINKTO is global, sequential is no
          CALL GXMAKE0(LINKTO(IDMN),NUM2,'LINKTO')
          CALL GXMAKE0(L0ADDQ(IDMN),NUM,'ADDQ')
        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 p
          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.GT.0.AND.CVAL.EQ.'GXO') THEN ! found ANGLED-OUT on this processo
            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.NE.' ') 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.EQ.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.EQ.'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.EQ.'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
          ENDIF
C
          CVAL=' '; CALL GETSDC(COBNAM,'INLET',CVAL)
          IF(CVAL.EQ.'GXI') THEN ! found ANGLED-IN
C... Deal with Angled inlets, and cartesian velocities in polar
            IF(II.GT.0) IIN=IIN+1
C... find parent blockage if set
            BLOCKNAM=' '; IBLK=0
            CALL GETSDC(COBNAM,'BLOCKAGE',BLOCKNAM) ! get name of parent
            IF(BLOCKNAM.NE.' ') 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.EQ.BLOCKNAM) THEN
                  IF(LPAR) THEN
                    IBLK=GD_INDPAT(III,1)
                  ELSE
                    IBLK=III
                  ENDIF
                  EXIT
                ENDIF
              ENDDO
            ENDIF
            IF(II.GT.0) F(L0BLKIN(IDMN)+IIN)=FLOAT(IBLK)
            BLOCKNAM=' '; IBLK=0
            CALL GETSDC(COBNAM,'LINKTO',BLOCKNAM) ! get name of linked
            IF(BLOCKNAM.NE.' ') THEN
              IF(BLOCKNAM.EQ.'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.EQ.'GXI') THEN
                    IBLK=III
                    EXIT
                  ENDIF
                ENDDO
                IF(IBLK.EQ.0) THEN
                  CALL WRIT40('No NEXT ANGLED-IN found for '
     1                                                     //OBJNAM(II))
                ENDIF
              ELSEIF(BLOCKNAM.EQ.'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.EQ.'GXI') THEN
                    IBLK=III
                    EXIT
                  ENDIF
                ENDDO
                IF(IBLK.EQ.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.EQ.BLOCKNAM) THEN
                    IBLK=III
                    EXIT
                  ENDIF
                ENDDO
                IF(IBLK.EQ.0) THEN
                  CALL WRIT40('Cannot find linked ANGLED-IN for '
     1                                                     //OBJNAM(II))
                ENDIF
              ENDIF
              IF(IBLK.EQ.0) THEN ! error finding linked ANGLED-IN
                CALL SET_ERR(574,'Cannot find linked ANGLED-IN for '
     1                                                   //OBJNAM(II),1)
              ENDIF
              IF(II.GT.0) THEN
                QADD=0.0
                CALL GETSDR(OBJNAM(II),'QADD',QADD) ! get additional heat
                F(L0ADDQ(IDMN)+IIN)=QADD
              ENDIF
            ENDIF
            IINL=ITWO(IP,IIN,LPAR)
            F(LINKTO(IDMN)+IINL)=FLOAT(IBLK)
C... set inlet density
            IF(RHO1.LE.0.0) THEN ! initialise inlet density
              RHOIN=1.0
            ELSE
              RHOIN=RHO1
            ENDIF
            CALL GETSDR(COBNAM,'DENIN',RHOIN)  ! get density at inlet
            IF(II.GT.0) F(L0DENIN(IDMN)+IIN)=RHOIN
C... Set inlet area ratio
            ARAT=1.0
            CALL GETSDR(COBNAM,'ARATI',ARAT)
            IF(II.GT.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.gt.0) THEN
C... for -ve total velocity, set density to GRND, to flag in-cell density
              IF(TOTVEL.LT.0.0) F(L0DENIN(IDMN)+IIN)=GRND
            ENDIF
            IF(II.GT.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.GT.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
            IF(QNE(VOLFLOW,-999.0).OR.QNE(MASFLOW,-999.0).OR.
     1                                      GRN(ENUT)) THEN ! need area
C... get patch limits, and sum surface area
              TOTAREA=0.0
              L0FACZ0=L0FACZ
              IF(II.LT.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(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(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
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
 1002         CONTINUE
              IF(LPAR) CALL GLSUM(TOTAREA) ! sum over all processors
              IF(II.GT.0) F(L0TOTAREA(IDMN)+IIN)=TOTAREA
C.... deduce total velocity from flowrate/area
              IF(QNE(VOLFLOW,-999.0).AND.II.GT.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.LT.0.0) F(L0DENIN(IDMN)+IIN)=GRND
              ELSEIF(QNE(MASFLOW,-999.0).AND.II.GT.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.GT.0) THEN
C... Inlet values for Turbulence
              TURBIN=-999.
              CALL GETSDR(COBNAM,'TURBIN',TURBIN)
              IF(QNE(TURBIN,-999.0)) THEN
                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)
              ENDIF
            ENDIF
            IF(.NOT.ONEPHS) THEN
C... set inlet density
              IF(RHO2.LE.0.0) THEN ! initialise inlet density
                RHOIN=1000.0
              ELSE
                RHOIN=RHO2
              ENDIF
              CALL GETSDR(COBNAM,'DENIN2',RHOIN)  ! get density at inlet
              IF(II.GT.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.GT.0) THEN
C... for -ve total velocity, set density to -1, to flag in-cell density
                IF(TOTVEL.LT.0.0) F(L0DENIN2(IDMN)+IIN)=-1.0
              ENDIF
              IF(II.GT.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.GT.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 r
C... get patch limits, and sum surface area
                TOTAREA=0.0
                IF(II.LT.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 stor
                  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.GT.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.LT.0.0) F(L0DENIN2(IDMN)+IIN)=-1.0
                ELSEIF(QNE(MASFLOW,-999.0).AND.II.GT.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(IDMN.LT.NUMDMN) THEN
          IDMN=IDMN+1
          GO TO 100
        ENDIF
C... reset indices for domain 1
        IF(NUMDMN.GT.1) THEN
          IDMN=1
          CALL INDXDM
        ENDIF
C... set index for SURN
        ISURN=LBNAME('SURN')
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.EQ.13.AND.ISC.EQ.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).GT.0
        IF(DOIT) THEN
          IPT=MAX(IOUT,IIN)
          DOIT=.FALSE.
          IF(INDVAR.EQ.P1.OR.INDVAR.EQ.P2.OR.INDVAR.EQ.R1.OR.
     1                                                INDVAR.EQ.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.EQ.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
            IBLK=NINT(F(L0BLKOUT(IDMN)+IPT))
            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.GT.0.0) THEN
                    IF(INDVAR.EQ.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.EQ.13.AND.ISC.EQ.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.GT.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(NCELL.GT.0) THEN
              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(NCELL.GT.0) THEN
              RELX=F(L0RL2(IDMN)+IPT)
              L0V=L0VL2(IDMN)+(IPT-1)*NCELL
            ENDIF
            LBVOUT=LBNAME('VOU2')
          ENDIF
          IF(LBVOUT.GT.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.EQ.U1.OR.INDVAR.EQ.U2) THEN  ! set direction vector
              CALL VECTOR(DIR,1.,0.,0.)
            ELSEIF(INDVAR.EQ.V1.OR.INDVAR.EQ.V2) THEN
              CALL VECTOR(DIR,0.,1.,0.)
            ELSEIF(INDVAR.EQ.W1.OR.INDVAR.EQ.W2) THEN
              CALL VECTOR(DIR,0.,0.,1.)
            ENDIF
            VEL=DOT(VELOUT,DIR) ! get velocity magnitude in velocity direction
          ELSEIF(QEQ(TOTVELO,-9999.)) THEN ! velocity is deduced
            L0MIN = L0PVAR(KPVMAS,IPAT,IZ)
            L0MINH= L0PVAR(KPVMAS,IPAT,IZ+1)
          ENDIF
          L0XU=L0F(XU2D)
          L0XG=L0F(XG2D)
          IBLK=NINT(F(L0BLKOUT(IDMN)+IPT))    ! get parent blockage patch number
          IDIR=(INDVAR-1)/2
          IPW=0
          IF(NCELL.GT.0) THEN
            CALL GETPAT(IPAT,IREG,RTYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
            IZADD=(IZ-JZF)*(JXL-JXF+1)*(JYL-JYF+1)
          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.GT.0.0) THEN
                  IF(INDVAR.GE.U1.AND.INDVAR.LE.W1) THEN ! momentum sources
C... VAL set to velocity component in relevant direction
                    IF(QEQ(TOTVELO,-999.)) THEN ! velocity vector has been set
                      IF(CARTES.OR.INDVAR.GT.V2) THEN
                        F(L0VAL+I)=VEL
                      ELSE  ! Polar U1-V2
                        UC=VELOUT(1)  ! set cartesian components
                        VC=VELOUT(2)
                        IF(INDVAR.LT.V1) THEN ! U1, U2
                          F(L0VAL+I)=UC*COS(F(L0XU+I))-VC*SIN(F(L0XU+I))
                        ELSE                  ! V1, V2
                          F(L0VAL+I)=UC*SIN(F(L0XG+I))+VC*COS(F(L0XG+I))
                        ENDIF
                      ENDIF
                    ELSEIF(QEQ(TOTVELO,-9999.)) THEN
C... deduce velocity from Vel = Mass_flow/(density*Area*Area_ratio)
                      L0MASS=ITWO(L0MIN,L0MINH,IPLUS.EQ.0)
                      VEL=F(L0MASS+IPW)/(ARAT*CAREA*F(L0DEN+I)+TINY)
                      IF(NCELL.GT.0) THEN
                        IV=IPW+IZADD
                        VEL=(1.-RELX)*F(L0V+IV)+RELX*VEL
                        F(L0V+IV)=VEL
                      ENDIF
                      IF(LBVOUT.GT.0) F(L0VOUT+I)=VEL
                      F(L0VAL+I)=NORML(IDIR)*VEL
                    ELSE ! velocity normal to surface has been set
                      F(L0VAL+I)=TOTVELO*NORML(IDIR) ! resolve into grid directi
                    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.GT.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.EQ.U1.OR.INDVAR.EQ.U2) THEN  ! set direction vector
              CALL VECTOR(DIR,1.,0.,0.)
            ELSEIF(INDVAR.EQ.V1.OR.INDVAR.EQ.V2) THEN
              CALL VECTOR(DIR,0.,1.,0.)
            ELSEIF(INDVAR.EQ.W1.OR.INDVAR.EQ.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
          IF(.NOT.CARTES) THEN
            L0XU=L0F(XU2D)
            L0XG=L0F(XG2D)
            IF(IURVAL.NE.0) L0RAD=L0F(RG2D)
          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(.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.EQ.W1.OR.INDVAR.EQ.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.GT.0.0) THEN ! found surface
                IF(INDVAR.EQ.P1.OR.INDVAR.EQ.P2.OR.INDVAR.EQ.R1.OR.
     1                                                INDVAR.EQ.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.GT.0.0) THEN ! fixed flow
                    IF(VELMASS.GE.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
                ELSEIF(INDVAR.EQ.ISURN) 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
                ELSE ! other sources
                  IF(INDVAR.GE.U1.AND.INDVAR.LE.W2) THEN ! momentum sources
C... VAL set to velocity component in relevant direction
                    IF(CARTES.OR.INDVAR.GT.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(INDVAR.LT.V1) THEN ! U1, U2
                        F(L0VAL+I)=UC*COS(F(L0XU+I))-VC*SIN(F(L0XU+I))
C... U is cartesian velocity component here, so negate Earth action for IURVAL n
                        IF(IURVAL.NE.0) THEN
                          F(L0VAL+I)=F(L0VAL+I)*F(L0RAD+IY)**IURVAL
                        ENDIF
                      ELSE                  ! V1, V2
                        F(L0VAL+I)=UC*SIN(F(L0XG+I))+VC*COS(F(L0XG+I))
                      ENDIF
                    ENDIF
                  ELSEIF(INDVAR.EQ.KE) THEN ! KE
                    F(L0VAL+I)=F(L0KEIN(IDMN)+IPT)
                  ELSEIF(INDVAR.EQ.EP.OR.INDVAR.EQ.LBOMEG.OR.INDVAR.
     1                            EQ.LBET) THEN ! EP, OMEG, ET
                    F(L0VAL+I)=F(L0EPIN(IDMN)+IPT)
                  ELSEIF(INDVAR.EQ.LBKP) THEN ! KP
                    F(L0VAL+I)=F(L0KEIN(IDMN)+IPT)/(1+.25)
                  ELSEIF(INDVAR.EQ.LBKT) THEN ! KT
                    GKTIN=0.25*F(L0KEIN(IDMN)+IPT)/(1+.25)
                  ENDIF
                ENDIF
              ENDIF
            ENDDO
          ENDDO
        ENDIF
      ELSEIF(IGR.EQ.13.AND.ISC.EQ.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.GT.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.EQ.U1.OR.INDVAR.EQ.U2) THEN  ! set direction vector
            CALL VECTOR(DIR,1.,0.,0.)
          ELSEIF(INDVAR.EQ.V1.OR.INDVAR.EQ.V2) THEN
            CALL VECTOR(DIR,0.,1.,0.)
          ELSEIF(INDVAR.EQ.W1.OR.INDVAR.EQ.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.NE.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.EQ.P1.OR.INDVAR.EQ.P2.OR.INDVAR.EQ.R1.OR.
     1                             INDVAR.EQ.R2.OR.INDVAR.EQ.ISURN) 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.NE.ISURN) THEN
                    IF(VELMASS.GT.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
                ELSE ! other sources
                  IF(INDVAR.GE.U1.AND.INDVAR.LE.W2) THEN ! momentum sources
C... VAL set to velocity component in relevant direction
                    IF(INDVAR.GT.V2) THEN ! Polar W1/W2
                      F(L0VAL+I)=VEL
                    ELSE  ! Polar U1-V2
                      IF(INDVAR.LT.V1) THEN ! U1, U2
                        F(L0VAL+I)=VELIN(1)*COS(F(L0XG+I))
     1                                          -VELIN(2)*SIN(F(L0XG+I))
C... U is cartesian velocity component here, so negate Earth action for
C    IURVAL ne. 0
                        IF(IURVAL.NE.0) THEN
                          F(L0VAL+I)=F(L0VAL+I)*F(L0RAD+IY)**IURVAL
                        ENDIF
                      ELSE                  ! V1, V2
                        F(L0VAL+I)=VELIN(1)*SIN(F(L0XG+I))
     1                                          +VELIN(2)*COS(F(L0XG+I))
                      ENDIF
                    ENDIF
                  ENDIF
                ENDIF
              ENDIF
            ENDDO
          ENDDO
        ENDIF
      ELSEIF(IGR.EQ.19.AND.ISC.EQ.7) THEN
C   * ------------------- SECTION 7 ---- Finish of sweep.
        IF(NGXIN(IDMN).EQ.0) RETURN ! nothing to do if no ANGLED-INs
C... Pick up sources for linked ANGLED-INs and make settings
        DO I=1,ILIM  ! loop looking for ANGLED-INs
          IF(LPAR) THEN
            COBNAM=GD_OBJNAM(I)
            CVAL=' '; CALL GETSDC(COBNAM,'INLET',CVAL)   ! try angled-in flag
            IF(CVAL.NE.'GXI') CYCLE ! not an ANGLED-IN
            IPT=I ! global patch number
            II=GD_INDPAT(IPT,1) ! get local patch number. -ve means not on this
          ELSE
            II=I
            IF(NINT(F(L0GXI(IDMN)+II)).EQ.0) CYCLE ! not an ANGLED-IN
            IPT=NINT(F(L0GXI(IDMN)+II)) ! sequence no of angled-in
          ENDIF
          IF(F(LINKTO(IDMN)+IPT).NE.0.) THEN ! this one is linked
            IR=NINT(F(LINKTO(IDMN)+IPT)) ! link flag
            IF(IR.GT.0) 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).EQ.'^') PATNM=GD_NAMPAT(ISOR)(2:)
              ISOR = INDPATCH(PATNM) ! local patch index from global name
            ENDIF
            IF(ISOR.LT.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.LT.0.0.AND.QNE(RHO1,GRND5).AND.ITEM1.EQ.0) THEN
              CALL SET_ERR(573,
     1          'Linked ANGLED-IN cannot use current density formula',1)
            ENDIF
            IF(.NOT.ONEPHS) THEN
              IF(ISOR.LT.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.LT.0.0.AND.QNE(RHO2,GRND5).AND.ITEM2.EQ.0) THEN
                CALL SET_ERR(573,
     1          'Linked ANGLED-IN cannot use current density formula',1)
              ENDIF
            ENDIF
            IF(ITAR.GT.0) ITR=NINT(F(L0GXI(IDMN)+ITAR)) ! sequence number of tar
            DO IPHI=14,NPHI ! loop over scalars
              IF(SOLVE(IPHI)) THEN
                IF(ISOR.LT.0)THEN
                  GSOR=0.0; GCO=0.0
                ELSE
                  CALL GETCV(ISOR,IPHI,GCO,GVAL) ! get C V from COVAL
                  IF(QEQ(GCO,-999.)) CYCLE ! no source was set so skip
                  CALL GETSO(ISOR,IPHI,GSOR) ! get source of variable in 'donor'
                ENDIF
                IF(LPAR) CALL GLSUM(GSOR)  ! sumover all processors
                IF(ITAR.LT.0) CYCLE ! target not on this processor
                IF(ONEPHS.OR.FLUID1(IPHI)) THEN
                  SMASS=SMASS1
                ELSE
                  SMASS=SMASS2
                ENDIF
                IF(NAME(IPHI).EQ.'TEM1') THEN ! Phase 1 temperature
                  IF(CP1.LT.0.0) CALL SET_ERR(572,
     1           'Linked ANGLED-IN cannot use variable specific heat',1)
                  QADD=F(L0ADDQ(IDMN)+ITR) ! additional heat
                  GAVE=(GSOR-QADD)/SMASS/CP1-TEMP0 ! subtract as SMASS is -ve!
                  T1AVE=GAVE
                ELSEIF(NAME(IPHI).EQ.'TEM2') THEN ! Phase 2 temperaure
                  IF(CP2.LT.0.0) CALL SET_ERR(572,
     1           'Linked ANGLED-IN cannot use variable specific heat',1)
                  QADD=F(L0ADDQ(IDMN)+ITR)
                  GAVE=(GSOR-QADD)/SMASS/CP2-TEMP0
                  T2AVE=GAVE
                ELSEIF(IPHI.EQ.14.OR.IPHI.EQ.15) THEN ! enthalpy
                  QADD=F(L0ADDQ(IDMN)+ITR)
                  GAVE=(GSOR-QADD)/SMASS
                ELSE    ! other scalars
                  GAVE=GSOR/SMASS ! mass-average scalar at 'donor'
                  IF(IPHI.EQ.C1) THEN
                    C1AVE=GAVE
                  ELSEIF(IPHI.EQ.C2) THEN
                    C2AVE=GAVE
                  ENDIF
                ENDIF
                CALL CORCOVAL(ITAR,IPHI,GCO,GAVE) ! set value back
              ENDIF
            ENDDO
            IF(ITAR.LT.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=-999.
              CALL GETSDR(OBJNAM(ITAR),'TURBIN',TURBIN)
              IF(QNE(TURBIN,-999.0)) THEN
                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
          ENDIF
        ENDDO
      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
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/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 '/phoeclos/d_includ/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)
      REAL NORML(3)
      LOGICAL DOIT,EQZ
C
      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 = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
      IF(IPBC.GT.0) THEN        ! dealing with cut cell
        IF(IBLK.EQ.0) THEN      ! no blockage specified - any will do
          DOIT=F(L0PAT+I).GT.0.0
        ELSE                    ! check for particular blockage
          DOIT=NINT(F(L0PAT+I)).EQ.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.LT.NX) ! neighbour indices - East
        INW=ITWO(I-NY, I,IX.GT. 1) !                   - West
        INN=ITWO(I+1,  I,IY.LT.NY) !                   - North
        INS=ITWO(I-1,  I,IY.GT. 1) !                   - South
        INH=ITWO(I+NXNY,I,IZ.LT.NZ) !                  - High
        INL=ITWO(I-NXNY,I,IZ.GT. 1) !                  - Low
        IF(F(L0PAT+INE).GT.0.0) THEN     ! solid to east
          IF(IBLK.EQ.0) THEN ! any blockage will do
            DOIT=.TRUE.
          ELSE               ! only next to blockage IBLK
            DOIT=NINT(F(L0PAT+INE)).EQ.IBLK
          ENDIF
          IF(DOIT) THEN
            IJKDME=IJKDM+NY              ! neighbour index
            IPBCE = IPBSEQ(IDMN,IJKDME)  ! neighbour cut-cell index
            IF(IPBCE.GT.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.GT.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).GT.0.0) THEN ! solid to west
          IF(IBLK.EQ.0) THEN
            DOIT=.TRUE.
          ELSE
            DOIT=NINT(F(L0PAT+INW)).EQ.IBLK
          ENDIF
          IF(DOIT) THEN
            IJKDMW=IJKDM-NY
            IPBCW = IPBSEQ(IDMN,IJKDMW)
            IF(IPBCW.GT.0) THEN
              ICUTIDXW = ISHPB(IDMN,IPBCW)
              AREA=F(PBC_DAKAR(3)+ICUTIDXW)
            ELSE
              AREA=F(I3DAEX+IJKDMW)
            ENDIF
            IF(AREA.GT.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).GT.0.0) THEN ! solid to north
          IF(IBLK.EQ.0) THEN
            DOIT=.TRUE.
          ELSE
            DOIT=NINT(F(L0PAT+INN)).EQ.IBLK
          ENDIF
          IF(DOIT) THEN
            IJKDMN=IJKDM+1
            IPBCN = IPBSEQ(IDMN,IJKDMN)
            IF(IPBCN.GT.0) THEN
              ICUTIDXN = ISHPB(IDMN,IPBCN)
              AREA=F(PBC_DAKAR(2)+ICUTIDXN)
            ELSE
              AREA=F(I3DANY+IJKDM)
            ENDIF
            IF(AREA.GT.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).GT.0.0) THEN ! solid to south
          IF(IBLK.EQ.0) THEN
            DOIT=.TRUE.
          ELSE
            DOIT=NINT(F(L0PAT+INS)).EQ.IBLK
          ENDIF
          IF(DOIT) THEN
            IJKDMS=IJKDM-1
            IPBCS = IPBSEQ(IDMN,IJKDMS)
            IF(IPBCS.GT.0) THEN
              ICUTIDXS = ISHPB(IDMN,IPBCS)
              AREA=F(PBC_DAKAR(1)+ICUTIDXS)
            ELSE
              AREA=F(I3DANY+IJKDMS)
            ENDIF
            IF(AREA.GT.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).GT.0.0) THEN ! solid to high
          IF(IBLK.EQ.0) THEN
            DOIT=.TRUE.
          ELSE
            DOIT=NINT(F(L0PAT+INH)).EQ.IBLK
          ENDIF
          IF(DOIT) THEN
            IJKDMH=IJKDM+NXNY
            IPBCH = IPBSEQ(IDMN,IJKDMH)
            IF(IPBCH.GT.0) THEN
              ICUTIDXH = ISHPB(IDMN,IPBCH)
              AREA=F(PBC_DAKAR(6)+ICUTIDXH)
            ELSE
              AREA=F(I3DAHZ+IJKDM)
            ENDIF
            IF(AREA.GT.CAREA) THEN
              CAREA=AREA
              CALL VECTOR(NORML,0.,0.,-1.)
            ENDIF
          ENDIF
        ENDIF
        IF(F(L0PAT+INL).GT.0.0) THEN ! solid to low
          IF(IBLK.EQ.0) THEN
            DOIT=.TRUE.
          ELSE
            DOIT=NINT(F(L0PAT+INL)).EQ.IBLK
          ENDIF
          IF(DOIT) THEN
            IJKDML=IJKDM-NXNY
            IPBCL = IPBSEQ(IDMN,IJKDML)
            IF(IPBCL.GT.0) THEN
              ICUTIDXL = ISHPB(IDMN,IPBCL)
              AREA=F(PBC_DAKAR(5)+ICUTIDXL)
            ELSE
              AREA=F(I3DAHZ+IJKDML)
            ENDIF
            IF(AREA.GT.CAREA) THEN
              CAREA=AREA
              CALL VECTOR(NORML,0.,0.,1.)
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      IF((INDVAR.EQ.W1.OR.INDVAR.EQ.W2).AND.EQZ(CAREA).AND.IZ.LT.NZ)THEN
C... If nothing found for W1, check next plane as well in H direction
        IJKDM=IJKDM+NXNY          ! increment 3D cell number
        IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
        L0PAT1=L0PAT+NXNY         ! increment patch-number store index
        IF(IPBC.GT.0) THEN        ! dealing with cut cell
          IF(IBLK.EQ.0) THEN
            DOIT=F(L0PAT1+I).GT.0.0
          ELSE
            DOIT=NINT(F(L0PAT1+I)).EQ.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(IZ.LT.NZ-1) THEN  ! dealing with un-cut cell
          INH=ITWO(I+NXNY,I,IZ.LT.NZ-1)
          IF(F(L0PAT1+INH).GT.0.0) THEN ! solid to high
            IF(IBLK.EQ.0) THEN
              DOIT=.TRUE.
            ELSE
              DOIT=NINT(F(L0PAT1+INH)).EQ.IBLK
            ENDIF
            IF(DOIT) THEN
              CALL VECTOR(NORML,0.,0.,-1.)
              IJKDMH=IJKDM+NXNY
              IPBCH = IPBSEQ(IDMN,IJKDMH)
              IF(IPBCH.GT.0) THEN
                ICUTIDXH = ISHPB(IDMN,IPBCH)
                CAREA=F(PBC_DAKAR(6)+ICUTIDXH)
              ELSE
                CAREA=F(I3DAHZ+IJKDMH)
              ENDIF
            ENDIF
          ENDIF
        ENDIF
        IF(CAREA.GT.0) IPLUS=1
      ENDIF
      END
c