c
SUBROUTINE CIRCLE(GBE,GBN,GBH,XC,YC,XLC,YLC,DIA,ICIR,NY,NX,
1 NS,ATOT)
C-----------------------------------------------------------------
C
C CONSIDER NOMINAL PLANE : X HORIZONTAL, Y VERTICAL
C
C GBE : BLOCKAGE IN X DIRECTION, COMMON, DIMENSIONED (Y,X)
C GBN : BLOCKAGE IN Y DIRECTION, COMMON, DIMENSIONED (Y,X)
C GBH : CELL BLOCKAGE (Z-DIRECTION), COMMON, DIMENSIONED (Y,X)
C XC : DISTANCES TO CELL FACES IN X DIRECTION, (1-D)
C YC : DISTANCES TO CELL FACES IN Y DIRECTION, (1-D)
C XLC : X DISTANCE TO CIRCLE CENTER
C YLC : Y DISTANCE TO CIRCLE CENTER
C DIA : CIRCLE DIAMETER
C ICIR : INTEGER DETERMINING GRID TYPE ( 1=POLAR,2=CARTESIAN )
C NS : INTEGER DETERMINING PARTITIONING OF GRID CELLS
C ATOT : AREA OF CIRCLE RETURNED BY SUBROUTINE
C
COMMON/CIRCOM/ICR,XL,YL,RAD
DIMENSION GBE(NY,NX),GBN(NY,NX),GBH(NY,NX)
DIMENSION XC(NX)
DIMENSION YC(NY)
DIMENSION XU(5),YV(5)
LOGICAL INXP,INXM,INYP,INYM
LOGICAL NEZ
INTEGER XDIR,YDIR
DATA XDIR,YDIR/1,2/
DATA TINY/1.E-10/
C
ICR=ICIR
XL=XLC
YL=YLC + TINY
ATOT=0.0
C
RAD=DIA/2.0
IF(DIA.LE.TINY) THEN
CALL PUT_LINE(' Zero or negative diameter in call to circle',
1 .true.)
CALL SET_ERR(441,
* 'Error. Zero or negative diameter in call to circle.',1)
C STOP
ENDIF
C.... Initialise porosity arrays
DO 2 JX= 1,NX
DO 2 JY= 1,NY
GBE(JY,JX) = 0.0
GBN(JY,JX) = 0.0
GBH(JY,JX) = 0.0
2 CONTINUE
NSY=MAX0(1,NS)
NSX=NSY
NSY=IFIX(YC(NY)*FLOAT(NSY)*10./(DIA*FLOAT(NY)))
NSX=IFIX(XC(NX)*FLOAT(NSX)*10./(DIA*FLOAT(NX)))
NSY = MAX0(NSY,1)
NSX = MAX0(NSX,1)
DNSYX=FLOAT(NSY*NSX)
IF(ICIR.EQ.1) THEN
CONX=CDELTA(YL,XDIR,INXP)
IF(.NOT.INXP) CONX=XC(NX)
CONY=RAD
ELSE
CONX=RAD
CONY=RAD
ENDIF
XMIN=AMAX1(0.0,XLC-CONX)
XMAX=AMIN1(XC(NX),XLC+CONX)
YMIN=AMAX1(0.0,YLC-CONY)
YMAX=AMIN1(YC(NY),YLC+CONY)
XP=0.0
DO 1000 IX=1,NX
XM=XP
XP=XC(IX)
XB=0.0
IF(XP.LT.XMIN) GO TO 1000
IF(XM.GT.XMAX) GO TO 1000
SX=(XP-XM)/FLOAT(NSX)
XP1=XM
DO 2000 ISX=1,NSX
XM1=XP1
XP1=XM1 + SX
YDELP=CDELTA(XP1,YDIR,INYP)
YDELM=CDELTA(XM1,YDIR,INYM)
IF(.NOT.INYP.AND..NOT.INYM) GO TO 2000
YP=TINY
YCONP=YLC
YCONM=YLC
IF(ICIR.EQ.1) THEN
YCONP=YCONP*COS(XLC-XP1)
YCONM=YCONM*COS(XLC-XM1)
ENDIF
DO 1001 IY=1,NY
YM=YP
YP=YC(IY)
YB=0.0
IF(YP.LT.YMIN) GO TO 1001
IF(YM.GT.YMAX) GO TO 1001
YP1=YM
SY=(YP-YM)/FLOAT(NSY)
IF(ICIR.EQ.1) SY=SY*(YP+YM)
DO 2001 ISY=1,NSY
XB=0.0
YB=0.0
YM1=YP1
YP1=YM1 + SY
IF(ICIR.EQ.1) YP1=SQRT(SY + YM1*YM1)
XDELP=CDELTA(YP1,XDIR,INXP)
XDELM=CDELTA(YM1,XDIR,INXM)
IF(.NOT.INXM.AND..NOT.INXP) GO TO 2001
XU(1)=XM1
XU(2)=XM1
XU(3)=XP1
XU(4)=XP1
XU(5)=XP1
YV(1)=YM1
YV(2)=YP1
YV(3)=YP1
YV(4)=YM1
YV(5)=YM1
C
C.... Calculate cell area
C
CALL GAREA(XU,YV,GTAREA,ICIR)
IF(NEZ(XLC).OR.ICIR.EQ.2) GO TO 999
IF(XDELP.LE.-.999999) GO TO 7000
IF(ABS(XDELP).GT.1.0) GO TO 91
IF(ABS(XDELM).GT.1.0) XDELM=-1.0
GO TO 6
91 CONTINUE
IF(ABS(XDELM).GT.1.001) GO TO 2001
IF(XDELM.LT.-0.99999) XDELM=-1.0
IF(ABS(XDELP).GT.1.0) XDELP=-1.0
GO TO 3
999 CONTINUE
IF(INXP) GO TO 6
IF(INXM) GO TO 3
GO TO 2001
3 IF(INYP) GO TO 100
IF(INYM) GO TO 100
GO TO 2001
6 IF(INXM) GO TO 10
IF(INYP) GO TO 400
IF(INYM) GO TO 400
GO TO 2001
10 IF(INYP) GO TO 12
IF(INYM) GO TO 700
GO TO 2001
12 IF(INYM) GO TO 800
GO TO 900
100 XSE=XLC+XDELM
XSW=XLC-XDELM
IF(XSE.LE.XP1) GO TO 110
IF(XSW.GE.XP1) GO TO 2001
IF(XSW.LT.XM1) GO TO 160
YEN=YCONP+YDELP
IF(YEN.GT.YP1) GO TO 955
GO TO 135
160 YEN=YCONP+YDELP
IF(YEN.GE.YP1) GO TO 145
GO TO 155
110 IF(XSE.LE.XM1) GO TO 2001
IF(XSW.GE.XM1) GO TO 125
GO TO 115
400 XNE=XLC+XDELP
XNW=XLC-XDELP
IF(XNE.LE.XP1) GO TO 410
IF(XNW.GE.XP1) GO TO 2001
IF(XNW.GE.XM1) GO TO 435
YWS=YCONM-YDELM
IF(YWS.LE.YM1) GO TO 445
GO TO 455
410 IF(XNE.LE.XM1) GO TO 2001
IF(XNW.GE.XM1) GO TO 425
YWS=YCONM-YDELM
IF(YWS.LE.YM1) GO TO 755
GO TO 415
700 YWN=YCONM+YDELM
YWS=YCONM-YDELM
IF(YWN.LE.YP1) GO TO 710
IF(YWS.GE.YP1) GO TO 2001
IF(YWS.GE.YM1) GO TO 415
760 XSE=XLC+XDELM
IF(XSE.GE.XP1) GO TO 745
GO TO 755
710 IF(YWN.LE.YM1) GO TO 2001
IF(YWS.GE.YM1) GO TO 725
GO TO 115
900 YEN=YCONP+YDELP
YES=YCONP-YDELP
IF(YEN.LE.YP1) GO TO 910
IF(YES.GE.YP1) GO TO 2001
IF(YES.GE.YM1) GO TO 435
960 XNW=XLC-XDELP
IF(XNW.LE.XM1) GO TO 945
GO TO 955
910 IF(YEN.LE.YM1) GO TO 2001
IF(YES.GE.YM1) GO TO 925
GO TO 135
800 YES=YCONP-YDELP
IF(YES.GE.YM1) GO TO 400
YWN=YCONM+YDELM
IF(YWN.LE.YP1) GO TO 100
XNE=XLC+XDELP
XNW=XLC-XDELP
IF(XNW.LE.XM1.AND.XNE.LE.XP1) GO TO 760
XSW=XLC-XDELM
IF(XSW.GE.XM1) GO TO 960
C.... Solid cell
7000 CONTINUE
POR=1.0
ATOT=ATOT + GTAREA
XB=XP1 - XM1
YB=YP1 - YM1
GO TO 9000
C.... Triangle, bottom LH corner
115 XU(3)=XLC+XDELM
XU(4)=XU(3)
XU(5)=XU(4)
YV(2)=YCONM+YDELM
YV(3)=YM1
GO TO 8000
C.... Semi circle
125 GO TO 2001
C.... Triangle, bottom RH corner
135 YEN=YCONP+YDELP
IF(YEN.GT.YP1) GO TO 955
XU(1)=XLC-XDELM
XU(2)=XP1
YV(2)=YCONP+YDELP
YV(3)=YM1
YB=YV(2) - YM1
GO TO 8000
C.... Pentangle, bottom RH corner
145 XU(3)=XLC-XDELP
YV(2)=YCONM+YDELM
YV(4)=YP1
XB=XP1 - XU(3)
YB=YP1 - YM1
GO TO 8000
C.... Rectangle, bottom
155 YV(2)=YCONM+YDELM
YV(3)=YCONP+YDELP
YB=YV(3) - YM1
GO TO 8000
C.... Triangle, top RH corner
435 XU(1)=XLC-XDELP
XU(2)=XU(1)
YV(1)=YP1
YV(4)=YCONP-YDELP
YV(5)=YV(4)
YB=YP1 - YV(4)
XB=XP1 - XU(1)
GO TO 8000
C.... Semi circle
425 GO TO 2001
C.... Triangle, top LH corner
415 XU(3)=XLC+XDELP
XU(4)=XU(3)
XU(5)=XU(4)
YV(1)=YCONM-YDELM
YV(4)=YP1
YV(5)=YP1
XB=XU(3) - XM1
GO TO 8000
C.... Pentangle, top LH corner
445 XU(5)=XLC+XDELM
YV(4)=YCONP-YDELP
YB=YP1 - YV(4)
XB=XP1 - XM1
GO TO 8000
C.... Top rentangle
455 YV(1)=YCONM-YDELM
YV(4)=YCONP-YDELP
YV(5)=YV(4)
YB=YP1 - YV(4)
XB=XP1 - XM1
GO TO 8000
C.... LH rectangle
755 XU(3)=XLC+XDELP
XU(5)=XLC+XDELM
XU(4)=XU(5)
XB=XU(3) - XM1
GO TO 8000
C.... Semi circle
725 GO TO 2001
C.... Pentangle bottom LH corner
745 XU(3)=XLC+XDELP
YV(4)=YCONP+YDELP
XB=XU(3) - XM1
YB=YV(4) - YM1
GO TO 8000
C.... RH rectangle
955 XU(1)=XLC-XDELM
XU(2)=XLC-XDELP
XB=XP1 - XU(2)
YB=YP1 - YM1
GO TO 8000
C.... Semi circle
925 GO TO 2001
C.... Pentangle, top RH corner
945 XU(5)=XLC-XDELM
YV(1)=YCONM-YDELM
YB=YP1 - YM1
XB=XP1 - XM1
C.... Calculate obstacle area
8000 CALL GAREA(XU,YV,GOAREA,ICIR)
ATOT=ATOT + GOAREA
POR=(GOAREA/GTAREA)
9000 GBH(IY,IX)=GBH(IY,IX) + POR/DNSYX
IF(ISY.EQ.NSY) GBN(IY,IX)=GBN(IY,IX) + XB
IF(ISX.EQ.NSX) GBE(IY,IX)=GBE(IY,IX) + YB
IF(ISY.EQ.NSY.AND.ISX.EQ.NSX)GBN(IY,IX)=GBN(IY,IX)/(XP-XM)
IF(ISX.EQ.NSX.AND.ISY.EQ.NSY)GBE(IY,IX)=GBE(IY,IX)/(YP-YM)
2001 CONTINUE
1001 CONTINUE
2000 CONTINUE
1000 CONTINUE
DO 9005 IX=1,NX
DO 9005 IY=1,NY
GBE(IY,IX) = AMAX1(0.0,GBE(IY,IX))
GBN(IY,IX) = AMAX1(0.0,GBN(IY,IX))
GBH(IY,IX) = AMAX1(0.0,GBH(IY,IX))
GBE(IY,IX) = AMIN1(1.0,GBE(IY,IX))
GBN(IY,IX) = AMIN1(1.0,GBN(IY,IX))
GBH(IY,IX) = AMIN1(1.0,GBH(IY,IX))
IXM=MAX0(1,IX-1)
IYM=MAX0(1,IY-1)
IF(GBE(IY,IX).LE.1.E-5.AND.GBE(IY,IXM).LE.1.E-5) THEN
IF(GBN(IY,IX).LE.1.E-5.OR.GBN(IYM,IX).LE.1.E-5) THEN
GBH(IY,IX )=0.0
GBE(IY,IX )=0.0
GBN(IY,IX )=0.0
GBE(IY,IXM)=0.0
GBN(IYM,IX)=0.0
ENDIF
ENDIF
IF(GBN(IY,IX).LE.1.E-5.AND.GBN(IYM,IX).LE.1.E-5) THEN
IF(GBE(IY,IX).LE.1.E-5.OR.GBE(IY,IXM).LE.1.E-5) THEN
GBH(IY,IX)=0.0
GBE(IY,IX)=0.0
GBN(IY,IX)=0.0
GBE(IY,IXM)=0.0
GBN(IYM,IX)=0.0
ENDIF
ENDIF
9005 CONTINUE
END
C****************************************************************
c
SUBROUTINE GAREA(GX,GY,GAR,JOPT)
C------------------------------------------------------------------
DIMENSION GX(5),GY(5)
CHARACTER BUFF*80
C
IF(JOPT.EQ.1) GO TO 10
IF(JOPT.EQ.2) GO TO 20
WRITE(BUFF,1)
CALL PUT_LINE(BUFF,.true.)
1 FORMAT(10X,23H(AREA) OPTION UNDEFINED/10X,23(1H*))
CALL SET_ERR(442, 'Error. See reslt file.',1)
C STOP
C------------------------------------------------------------------
CHAPTER 1 POLAR CO-ORDINATES
C------------------------------------------------------------------
10 CONTINUE
GA1=GY(1)*GY(4)*SIN(GX(4)-GX(1))
GA2=GY(1)*GY(5)*SIN(GX(5)-GX(1))
GA3=GY(5)*GY(4)*SIN(GX(4)-GX(5))
GA4=GY(2)*GY(3)*SIN(GX(3)-GX(2))
GA5=GY(1)*GY(3)*SIN(GX(3)-GX(1))
GA6=GY(1)*GY(2)*SIN(GX(2)-GX(1))
GA7=GY(3)*GY(4)*SIN(GX(4)-GX(3))
GAR=0.5*(ABS(GA1-GA2-GA3)+ABS(GA4-GA5+GA6)+ABS(GA1-GA5-GA7))
RETURN
C------------------------------------------------------------------
CHAPTER 2 CARTESIAN CO-ORDINATES
C------------------------------------------------------------------
20 CONTINUE
GB1=GY(1)*(GX(3)-GX(2))+GY(2)*(GX(1)-GX(3))+GY(3)*(GX(2)-GX(1))
GB2=GY(1)*(GX(3)-GX(4))+GY(3)*(GX(4)-GX(1))+GY(4)*(GX(1)-GX(3))
GB3=GY(1)*(GX(4)-GX(5))+GY(4)*(GX(5)-GX(1))+GY(5)*(GX(1)-GX(4))
GAR=0.5*(ABS(GB1)+ABS(GB2)+ABS(GB3))
END
C****************************************************************
FUNCTION CDELTA(Z,IDIR,INSIDE)
C---------------------------------------------------------------------
COMMON/CIRCOM/ICIR,XLC,YLC,RAD
LOGICAL INSIDE
INSIDE=.TRUE.
IF(ICIR.EQ.2) THEN
C.... Calculate X-direction intercepts (ie X=XLC +or- CDELTA)
IF(IDIR.EQ.1) THEN
CDELTA=-Z*Z-YLC*YLC+RAD*RAD+2.*YLC*Z
IF(CDELTA.LT.0.0) THEN
INSIDE=.FALSE.
CDELTA=0.0
RETURN
ENDIF
CDELTA=SQRT(CDELTA)
C.... Calculate Y-direction intercepts (ie Y=YLC +or- CDELTA)
ELSEIF(IDIR.EQ.2) THEN
CDELTA=-Z*Z-XLC*XLC+RAD*RAD+2.*XLC*Z
IF(CDELTA.LT.0.0) THEN
INSIDE=.FALSE.
CDELTA=0.0
RETURN
ENDIF
CDELTA=SQRT(CDELTA)
ENDIF
C.... Polar co-ordinates
ELSEIF(ICIR.EQ.1) THEN
C.... Calculate X-direction intercepts (ie X=XLC +or- CDELTA)
IF(IDIR.EQ.1) THEN
CDELTA=(Z*Z+YLC*YLC-RAD*RAD)/(2.*YLC*Z)
IF(ABS(CDELTA).GT.1.0) THEN
INSIDE=.FALSE.
CDELTA=AMIN1(-1.0,CDELTA)
RETURN
ENDIF
CDELTA=ACOS(CDELTA)
C.... Calculate Y-direction intercepts (ie Y=YLC +or- CDELTA)
ELSEIF(IDIR.EQ.2) THEN
CDELTA=(YLC*COS(XLC-Z))**2-YLC*YLC+RAD*RAD
IF(CDELTA.LT.0.0) THEN
INSIDE=.FALSE.
CDELTA=0.0
RETURN
ENDIF
CDELTA=SQRT(CDELTA)
ENDIF
ENDIF
END
C***********************************************************************
c
SUBROUTINE SLOPE(PC,PN,PH,YVDIST,ZWDIST,IBCOPT,THETA,CARTES
1 ,YORIG,ZORIG,IYF,IYL,IZF,IZL,NY,NZ,SURFA)
C-----------------------------------------------------------------------
C SET BLOCKAGE OPTION (IBCOPT) TO :-
C 1, FOR POSITIVE GRADIENT OF SLOPE : SOUTH-HIGH BLOCKED
C 2, FOR POSITIVE GRADIENT OF SLOPE : NORTH-LOW BLOCKED
C 3, FOR NEGATIVE GRADIENT OF SLOPE : NORTH-HIGH BLOCKED
C 4, FOR NEGATIVE GRADIENT OF SLOPE : SOUTH-LOW BLOCKED
C THETA IS ANGLE FROM Z(OR X)-AXIS : 0
SUBROUTINE CNESRF(IOPT,YM1,YP1,ZM1,ZP1,YO,ZO,TAN,CAREA)
C-----------------------------------------------------------------------
C This subroutine is called by SUBROUTINE SLOPE.
C-----------------------------------------------------------------------
C.... Define functions
C
YLOC(DUM)= YO+(DUM-ZO)*TAN
ZLOC(DUM)= ZO+(DUM-YO)/TAN
C
C.... Calculate intercept locations
C
YLOC1 = YLOC( ZM1 )
YLOC2 = YLOC( ZP1 )
ZLOC1 = ZLOC( YM1 )
ZLOC2 = ZLOC( YP1 )
IF(IOPT.GT.2) THEN
YN = AMIN1( YP1 , YLOC1 )
YS = AMAX1( YM1 , YLOC2 )
ZH = AMIN1( ZP1 , ZLOC1 )
ZL = AMAX1( ZM1 , ZLOC2 )
ELSE
YN = AMIN1( YP1 , YLOC2 )
YS = AMAX1( YM1 , YLOC1 )
ZH = AMIN1( ZP1 , ZLOC2 )
ZL = AMAX1( ZM1 , ZLOC1 )
ENDIF
C
C.... Calculate surface area
C
DZ= ZH - ZL
DY= YN - YS
SLANTH= SQRT(DY*DY + DZ*DZ)
CAREA= (YN + YS) * SLANTH * 0.5
C
END
C***********************************************************************
c
SUBROUTINE SAREA(X,Y,AR)
C***********************************************************************
DIMENSION X(5),Y(5)
B1=Y(1)*(X(3)-X(2))+Y(2)*(X(1)-X(3))+Y(3)*(X(2)-X(1))
B2=Y(1)*(X(3)-X(4))+Y(3)*(X(4)-X(1))+Y(4)*(X(1)-X(3))
B3=Y(1)*(X(4)-X(5))+Y(4)*(X(5)-X(1))+Y(5)*(X(1)-X(4))
AR=0.5*(ABS(B1)+ABS(B2)+ABS(B3))
END
c