C ALGORITHM 625 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 4, C DEC., 1984, P. 453. SUBROUTINE REGION(XGRID, YGRID, NGDIMX, NGDIMY, BRANGE, NPDIM, REG 10 * BCOORD, SCLOCK, SARC, SLEVEL, GTYPE, XBOUND, YBOUND, PIECE, * BPTYPE, BNEIGH, BGRID, BPARAM, NBDIM, NBPTS, FAIL) C C DOMAIN PROCESSOR. WRITTEN BY JOHN R. RICE, PURDUE UNIVERSITY. C C THE PURPOSE OF THIS ALGORITHM IS TO RELATE A GENERAL 2 DIMENSIONAL C DOMAIN TO A RECTANGULAR GRID LAID OVER IT. THE DOMAIN BOUNDARY IS C GIVEN IN PARAMETRIC FORM AND THE PROGRAM PRODUCES ARRAYS WHICH C DESCRIBE THE RELATIONSHIP BETWEEN THE TWO OBJECTS. THESE ARRAYS C CONTAIN INFORMATION POINTING IN BOTH DIRECTIONS (FROM THE DOMAIN C GRID TO THE BOUNDARY DATA STRUCTURE AND VICE VERSA), THIS SHOULD C FACILITATE APPLICATIONS( E.G. NUMERICAL QUADRATURE, SURFACE C FITTING, DISCRETIZING PDE'S ) WHICH INVOLVE SUCH DOMAINS. C C NOTE THAT SOME VARIABLES HAVE DUPLICATE NAMES SO THEY CAN BE BOTH C IN COMMON AND IN THE ARGUMENT LIST OF REGION C C ------------------- INPUT INFORMATION FOR DOMAIN PROCESSOR -------- C C 1** OUTPUT LEVEL CONTROL ******* C SLEVEL = LEVEL = OUTPUT CONTROL SETTING C DETAILS GIVEN IN SUBROUTINE DOMAIN C C 2** GRID DEFINITION ************ C C NGDIMX,NGDIMY = NUMBER OF GRID LINES IN X AND Y COORDINATES C = ACTUAL DECLARED DIMENSIONS OF GTYPE ARRAY C NGRIDX = NGDIMX, NGRIDY = NGDIMY C XGRID(IX),YGRID(JY) FOR IX = 1 TO NGRIDX, JY = 1 TO NGRIDY C = ARRAYS CONTAINING THE GRID LINE COORDINATES C *** THESE MUST BE IN ASCENDING ORDER *** C C 3** BOUNDARY DEFINITION ******** C C NPDIM = NBOUND = NUMBER OF BOUNDARY PIECES C NBDIM = ARRAY DIMENSION FOR BOUNDARY POINTS C BCOORD = PARAMETERIZED DEFINITION OF THE BOUNDARY. C BCOORD(P,X,Y,IPIECE) GIVES THE X,Y VALUES OF THE C POINT ON PIECE IPIECE WITH PARAMETER VALUE = P. C SEE NOTE BELOW FOR DISCRETE BOUNDARY DEFINITION C C BRANGE(1,I),BRANGE(2,I) C = FIRST AND LAST VALUES OF PARAMETERS DEFINING C THE I-TH BOUNDARY PIECE C SCLOCK = CLOCK = SWITCH TO SPECIFY BOUNDARY ORIENTATION C .TRUE. MEANS BOUNDARY IS CLOCKWISE C .FALSE. MEANS BOUNDARY IS COUNTER-CLOCKWISE C SARC = ARC = .TRUE. MEANS DOMAIN IS AN ARC WITH NO INTERIOR C C --------------------------------------------------------------------- C --------- OUTPUT INFORMATION OF THE DOMAIN PROCESSOR IS IN TWO PARTS C ( THERE IS ALSO A FAILURE FLAG = FAIL ) C C PART 1 ********** GRID POINT TYPES ************ C C GTYPE(IX,JY) FOR IX = 1 TO NGRIDX, JY = 1 TO NGRIDY C THE VALUES IN THIS ARRAY GIVE THE TYPE OF THE GRID POINTS AND C INFORMATION ABOUT THEIR RELATION TO THE BOUNDARY. DETAILS AND C EXAMPLES ARE GIVEN IN THE SUBROUTINE DOMAIN. C POSSIBLE VALUE ARE( ASSUMING PACKING FACTOR IPACKB = 1000) C C = INTEGER OVER 1000 C GRID POINT IS NEXT TO THE BOUNDARY AND THE GTYPE VALUE IS C GTYPE = K + 1000*J WHERE C K IS THE INDEX OF THE LOWEST NUMBERED BOUNDARY NEIGHBOR C J = FOUR BITS TO NOTE LOCATION OF BOUNDARY POINTS C 0001 - BOUNDARY NEIGHBOR TO NORTH (NOON) C 0010 - BOUNDARY NEIGHBOR TO EAST (3 O'CLOCK) C 0100 - BOUNDARY NEIGHBOR TO SOUTH (6 O'CLOCK) C 1000 - BOUNDARY NEIGHBOR TO WEST (9 O'CLOCK) C IF GTYPE = 6011 THEN THIS POINT HAS BOUNDARY NEIGHBORS C TO THE SOUTH AND EAST, 11 IS LOWEST NUMBERED OF THEM C C = 999 C GRID POINT IS INTERIOR TO REGION AND NOT CLOSE TO BOUNDARY C C = INTEGER LESS THAN 999 C GRID POINT IS ALSO BOUNDARY PT., GTYPE IS ITS INDEX C C = 0 C GRID POINT IS EXTERIOR FAR FROM THE BOUNDARY C C = NEGATIVE INTEGER C GRID POINT IS EXTERIOR NEXT TO THE BOUNDARY, ITS LOCATION C RELATIVE TO THE BOUNDARY IS ENCODED AS FOR INTERIOR POINTS C NEAR THE BOUNDARY C C C PART 2 ********* BOUNDARY/GRID INTERSECTIONS ********** C C NBPTS = NBNDPT = NUMBER OF BOUNDARY POINTS ACTUALLY FOUND C BOUNDARY POINT INDEX RANGE = 1 TO NBNDPT+1 C THE FIRST BOUNDARY POINT IS ADDED TO THE LIST TO C MAKE IT CIRCULAR, HENCE THE +1 IN THE INDEX RANGE C XBOUND(I),YBOUND(I) = COORDINATES OF I-TH BOUNDARY POINT C PIECE(I) = INDEX OF BOUNDARY PIECE TO WHICH PT. BELONGS C SMALLEST NUMBER FOR CORNER POINTS C BPTYPE(I) = TYPE OF BOUNDARY POINT C = HORZ IF POINT IN ON A Y GRID LINE C = VERT X C = BOTH IF POINT IS ALSO A GRID POINT C = INTE IF POINT IS NOT ON A GRID LINE C HAPPENS ONLY FOR CORNERS OF THE BOUNDARY C = JUMP IF POINT PRECEEDS A HOLE C BNEIGH(I) = POINTER TO THE INTERIOR POINTS FROM THE I-TH C BOUNDARY POINT. SAME SCHEME IS USED TO ENCODE C DIRECTIONS AS FOR THE J PART OF GTYPE ABOVE C BGRID(I) = IX + IPACKB*JY WHEN PT. I IS IN GRID SQUARE IX,JY C BPARAM(I) = PARAMETER VALUE OF I-TH BOUNDARY POINT C C ********************************************************************* C ************** DOMAIN PROCESSOR ALGORITHM STRUCTURE *************** C C THERE ARE 2 BASIC PARTS TO THE ALGORITHM, PROCESSING THE BOUNDARY C AND TYPING THE GRID POINTS. AN OUTLINE FOLLOWS C C MARK ALL GRID POINTS AS EXTERIOR C C******* PART 1 *** PROCESS BOUNDARY C C LOCATE FIRST BOUNDARY POINT C C DO LOOP OVER BOUNDARY PIECES C 1 DO WHILE NOT AT END OF PIECE C 1 1 LOCATE NEXT INTERSECTION WITH GRID LINE( SUBROUTINE BWALK ) C 1 1 DETERMINE TYPE FOR INTERSECTION POINT FOUND C 1 1 CHECK FOR CHANGING BOUNDARY PIECE(SUBROUTINE CHANGE DOES IT) C 1 1 CHECK CONTINUITY OF BOUNDARY C 1 +--- C +------ C C CHECK CLOSING OF BOUNDARY C FINISH UP FIRST POINT C C******* PART 2 *** TYPE GRID POINTS AND MARK INTERIOR POINTS C C PASS OVER GRID MARKING THE GRID POINT TYPES(SUBROUTINE GVALUS ) C MARK INTERIOR POINTS( SUBROUTINE FILL ) C LOCATE FIRST INTERIOR POINT C MARK THE INTERIOR POINTS C SET POINTERS FROM BOUNDARY TO INTERIOR( SUBROUTINE NEIGH ) C C C **** NOTE. IF THE BOUNDARY IS DESCRIBED DISCRETELY BY THE SET OF C BOUNDARY-GRID INTERSECTIONS, THEN THE FOLLOWING FUNCTION PROGRAM C ALLOWS THE DOMAIN PROCESSOR TO BE APPLIED. IT ESSENTIALLY PROVIDES C LINEAR INTERPOLATION BETWEEN THE DISCRETE POINTS C C * FUNCTION BCOORD(P,X,Y,IPIECE) C * REAL XBOUND(***),YBOUND(***) C * LOGICAL STARTD C * DATA STARTD/.FALSE./ C * IF( STARTD ) GO TO 20 C * OBTAIN DATA FOR DISCRETE POINTS THROUGH A COMMON BLOCK C * OR BY READING FROM A FILE. THE READ CASE IS SHOWN. C * READ(INFILE,10) NBNDPT,(XBOUND(I),YBOUND(I),I=1,NBNDPT-1) C * 10 FORMAT( ****** ) C * XBOUND(NBNDPT) = XBOUND(1) C * YBOUND(NBNDPT) = YBOUND(1) C * C * OBTAIN COORDINATES FROM ARRAYS OF DISCRETE DATA C * 20 IP = P C * X = (P-IP)*XBOUND(IP+1) + (IP+1-P)*XBOUND(IP) C * Y = (P-IP)*YBOUND(IP+1) + (IP+1-P)*YBOUND(IP) C * RETURN C * END C C ***************************************************************** C ***** REGION IS AN INTERFACE SUBROUTINE FOR THE DOMAIN PROCESSING C MAIN SUBPROGRAM = DOMAIN. C REGION MUST BE THE FIRST ROUTINE CALLED FOR PROCESSING C MULTIPLE REGIONS( E.G. WITH HOLES) C THE FUNCTIONS OF THIS INTERFACE PROGRAM ARE AS FOLLOWS. C C 1. INITIALIZE OR COMPUTE VARIOUS CONSTANTS AND PLACE THEM IN C COMMON BLOCKS. THESE ARE C A. I/O CONSTANT = MOUTPT C B. NUMERICAL CONTROL CONSTANTS = EPSGRD,EPSTAN C C. 5 CHARACTER STRING CONSTANTS C D. COUNTERS = N1BND,N1BDPT = 1 C E. INTERIOR MARKER = NSIDE C F. PACKING FACTOR = IPACKB( SHOULD BE NSIDE+1 ) C C 2. CHECK THAT THE GRID IS PROPERLY DEFINED. C C 3. MOVE VARIABLES FROM SUBROUTINE ARGUMENT LIST INTO COMMON. C A. DIMENSION SPECIFICATIONS = NGDIMX,NGDIMY,NBOUND C B. CONTROL VARIABLES = LEVEL, CLOCKW, ARC C C 4. CALL DOMAIN. C C 5. SET FAILURE FLAG FAIL, FINAL NUMBER NBPTS OF BOUNDARY POINTS C C REAL XGRID(NGDIMX), YGRID(NGDIMY), BRANGE(2,NPDIM) LOGICAL SCLOCK, SARC, FAIL INTEGER SLEVEL DIMENSION GTYPE(NGDIMX,NGDIMY), XBOUND(NBDIM), YBOUND(NBDIM), * BPTYPE(NBDIM), BNEIGH(NBDIM), BGRID(NBDIM), BPARAM(NBDIM), * PIECE(NBDIM) C C ******** IF HOLE IS ALSO USED THEN THE CALLING PROGRAM ***** C ******** MUST HAVE THE FOLLOWING COMMON BLOCKS IN ORDER ***** C ******** TO CONFORM TO FORTRAN STANDARDS ***** C EXTERNAL BCOORD INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C C CONSTANT CHARACTER STRINGS FOR NICE PRINT OUT INTEGER IQHORZ, IQVERT, IQBOTH, IQINTE, IQJUMP DATA IQHORZ, IQVERT, IQBOTH, IQINTE, IQJUMP /4HHORZ,4HVERT,4HBOTH, * 4HINTE,4HJUMP/ C C INITIALIZE COMMON BLOCKS FROM SUBROUTINE ARGUMENTS NGRIDX = NGDIMX NGRIDY = NGDIMY NBOUND = NPDIM LEVEL = SLEVEL CLOCKW = SCLOCK ARC = SARC C C INITIALIZE CONSTANT CHARACTER STRINGS FROM DATA STATEMENTS HORZ = IQHORZ VERT = IQVERT BOTH = IQBOTH INTER = IQINTE JUMP = IQJUMP C SET UNIT NUMBER FOR OUTPUT MOUTPT = 6 C PRINT GREETINGS C IF (LEVEL.GT.0) WRITE (MOUTPT,99999) 99999 FORMAT (///1H , 19(1H-)/1H , 18H DOMAIN PROCESSOR /1H , 19(1H-)// * 5X, 33H D O M A I N P R O C E S S O R ) C C SET MARKER FOR INTERIOR POINTS OF REGION C MUST EXCEED MAX NUMBER OF BOUNDARY POINTS NSIDE = 999 IPACKB = NSIDE + 1 C SET LIMIT ON THE QUEUE LENGTH IN ROUTINE FILL C QLIMIT = 250 SHOULD ALLOW WELL OVER 10,000 C INTERIOR POINTS IN THE DOMAIN QLIMIT = 250 C MARK THAT WE ARE NOT IN HOLE OR ARC ARC = .FALSE. INHOLE = .FALSE. C INITIALIZE BOUNDARY LIST COUNTERS N1BND = 1 N1BDPT = 1 C C SET GEOMETRY TOLERANCE PARAMETER EPSGRD C **************************************************************** C ****** ****** C ****** THIS IS A MACHINE AND PROBLEM DEPENDENT CONSTANT ****** C ****** ****** C **************************************************************** C C EPSGRD SHOULD BE LARGE ENOUGH TO INSULATE THE COMPUTATIONS FROM C MACHINE ROUND-OFF. ALL POINTS AND LINES WITHIN EPSGRD OF ONE C ANOTHER ARE ASSUMED TO BE EQUAL. IT IS PROBABLY SAFE TO TAKE C EPSGRD AS 50 UNITS IN THE LAST PLACE, IT MUST BE AT LEAST 20 C UNITS IN THE LAST PLACE. THE CONVERGENCE TEST IN SECANT IS C .2*EPSGRD C C ****** THIS PARAMETER IS RELATIVE TO THE DOMAIN SIZE ****** C EPSGRD = FACTOR*( DOMAIN WIDTH ) C C EPSGRD SHOULD BE SMALL ENOUGH SO THAT THE ACCURACY IN THE PROBLEM C SOLUTION IS NOT AFFECTED BY AN UNCERTAINTY IN THE GEOMETRY OF C EPSGRD. C C FACTOR = 1.E-8 IS APPROPRIATE FOR MOST LONG WORD LENGTH C MACHINES AND MOST PROBLEMS C 2.E-5 IS APPROPRIATE( BUT NOT FAIL-SAFE ) FOR C 32 BIT MACHINES. XWIDTH = XGRID(NGRIDX) - XGRID(1) YWIDTH = YGRID(NGRIDY) - YGRID(1) EPSGRD = 2.E-5*AMAX1(XWIDTH,YWIDTH) C C CHECK THAT THE GRIDS ARE PROPERLY DEFINED IF (NGRIDX.GE.2) GO TO 20 C FATAL ERROR, GRID TOO SMALL 10 CONTINUE C IF (LEVEL.GE.0) WRITE (MOUTPT,99998) 99998 FORMAT (/5(3H **), 39H FATAL ERROR, MUST HAVE AT LEAST 2 GRID, * 18H LINES IN X AND Y ) FATAL = .TRUE. RETURN 20 IF (NGRIDY.LE.1) GO TO 10 C C CHECK THAT THE GRID INCREASES IN BOTH DIRECTIONS C FIND THE MINIMUM GRID WIDTHS XGMIN = XWIDTH YGMIN = YWIDTH DO 30 I=2,NGRIDX XGMIN = AMIN1(XGMIN,XGRID(I)-XGRID(I-1)) 30 CONTINUE DO 40 I=2,NGRIDY YGMIN = AMIN1(YGMIN,YGRID(I)-YGRID(I-1)) 40 CONTINUE IF (AMIN1(XGMIN,YGMIN).GT.EPSGRD) GO TO 50 C C FATAL ERROR, HAVE ZERO OR NEGATIVE GRID WIDTH C IF (LEVEL.GE.0) WRITE (MOUTPT,99997) 99997 FORMAT (/5(3H **), 32H FATAL ERROR, X AND Y GRID LINES/9X, * 19H MUST BE INCREASING) FATAL = .TRUE. RETURN C C SET PARAMETER FOR TANGENCY TEST C SEE SUBROUTINE CHKTAN FOR DETAILS 50 EPSTAN = .1E0*AMIN1(XGMIN,YGMIN,15.E0*SQRT(EPSGRD)) C C CALL DOMAIN PROCESSOR CALL DOMAIN(XGRID, YGRID, NGDIMX, NGDIMY, BRANGE, NPDIM, BCOORD, * GTYPE, XBOUND, YBOUND, PIECE, BPTYPE, BNEIGH, BGRID, BPARAM, * NBDIM) C C SET FAILURE FLAG AND FINAL NUMBER OF BOUNDARY POINTS FAIL = FATAL C THESE NUMBERS WILL BE USED BY HOLE IF IT IS CALLED NBPTS = NBNDPT N1BND = NBOUND + 1 C CHECK FOR INTERIOR POINTS ON BOUNDARY OF GTYPE ARRAY C IF FOUND IT MEANS THERE IS SOMETHING WRONG DO 60 I=1,NGDIMX IF (GTYPE(I,1).EQ.NSIDE .OR. GTYPE(I,NGDIMY).EQ.NSIDE) GO TO 80 60 CONTINUE DO 70 J=1,NGDIMY IF (GTYPE(1,J).EQ.NSIDE .OR. GTYPE(NGDIMX,J).EQ.NSIDE) GO TO 80 70 CONTINUE C NO PROBLEM DETECTED RETURN C SERIOUS PROBLEM DETECTED 80 CONTINUE C IF (LEVEL.GE.0) WRITE (MOUTPT,99996) 99996 FORMAT (/5(3H **), 33H FATAL ERROR, FOUND INTERIOR PTS /4X, * 43HON THE EDGE OF THE GRID. PROBABLE CAUSES /4X, 10H BOUNDAR, * 30HY ORIENTATION IS WRONG /4X, 23H BOUNDARY STARTS WHER, * 24HE INTERIOR IS TOO THIN /4X, 29H BOUNDARY OSCILLATES TOO RA, * 17HPIDLY SOME PLACE ) FATAL = .TRUE. RETURN END SUBROUTINE DOMAIN(XGRID, YGRID, NGDIMX, NGDIMY, BRANGE, NPDIM, DOM 10 * BCOORD, GTYPE, XBOUND, YBOUND, PIECE, BPTYPE, BNEIGH, BGRID, * BPARAM, NBDIM) C C ***** THIS SUBROUTINE PROCESSES THE RECTANGULAR GRID AND SPECIFIED C BOUNDARY. IT APPLIES TO ONE CLOSED LOOP OR ARC OF THE BOUNDARY AND C MAY BE CALLED SEVERAL TIMES FOR A COMPLEX DOMAIN C C ------------------- INPUT INFORMATION FOR DOMAIN ------------------ C THE INPUT INFORMATION IS IN THE COMMON BLOCKS AND ARGUMENTS C SEE THE MAIN DRIVER REGION FOR MORE DETAILS C C 1** OUTPUT LEVEL CONTROL ******* C LEVEL = CONTROL SETTING, DETAILS GIVEN BELOW C C 2** GRID DEFINITION ************ C C NGDIMX,NGDIMY = NGRIDX,NGRIDY = GRID LINES IN X AND Y COORDINATES C XGRID(IX),YGRID(JY) FOR IX = 1 TO NGRIDX, JY = 1 TO NGRIDY C C 3** BOUNDARY DEFINITION ******** C C N1BND = NUMBER OF THE FIRST BOUNDARY PIECE C WILL DIFFER FROM 1 AFTER FIRST CALL C NPDIM = ARRAY DIMENSION FOR BOUNDARY PIECES C N1BDPT = NUMBER OF THE FIRST BOUNDARY POINT C WILL DIFFER FROM 1 AFTER FIRST CALL C NBDIM = ARRAY DIMENSION FOR BOUNDARY POINTS C BCOORD = PARAMETERIZED DEFINITION OF THE BOUNDARY. C BCOORD(P,X,Y,IPIECE) GIVES THE X,Y VALUES OF THE C POINT ON PIECE IPIECE WITH PARAMETER VALUE = P. C BRANGE(2,I) = FIRST AND LAST VALUES OF PARAMETERS DEFINING C THE I-TH BOUNDARY PIECE C CLOCKW = SWITCH TO SPECIFY BOUNDARY ORIENTATION C .TRUE. MEANS BOUNDARY IS CLOCKWISE C .FALSE.MEANS BOUNDARY IS COUNTER-CLOCKWISE C ARC = .TRUE. MEANS DOMAIN IS AN ARC WITH NO INTERIOR C C ----------- THE OUTPUT OF DOMAIN IS IN TWO PARTS ------------- C C 1** GRID SPECIFICATION ********** C C GTYPE(IX,JY) FOR IX = 1 TO NGRIDX, JY = 1 TO NGRIDY C THE VALUES IN THIS ARRAY GIVE THE TYPE OF THE GRID POINTS AND C INFORMATION ABOUT THEIR RELATION TO THE BOUNDARY. C THERE IS A PACKING FACTOR IPACKB WHICH IS NORMALLY 1000. FOR C VERY LARGE PROBLEMS, IPACKB AND RELATED CONSTANT NSIDE = C IPACKB - 1 MUST BE INCREASED SO THAT NSIDE .GT. NBNDPT. C POSSIBLE VALUE ARE( ASSUMING IPACKB = 1000) C C = INTEGER OVER 1000 C GRID POINT IS NEXT TO THE BOUNDARY AND THE GTYPE VALUE IS C GTYPE = K + 1000*J WHERE C K IS THE INDEX OF THE LOWEST NUMBERED BOUNDARY NEIGHBOR C MUST DOUBLE CHECK USE WHEN K = 1 C J = FOUR BITS TO NOTE LOCATION OF BOUNDARY POINTS C 0001 - BOUNDARY NEIGHBOR TO NORTH (NOON) C 0010 - BOUNDARY NEIGHBOR TO EAST (3 O'CLOCK) C 0100 - BOUNDARY NEIGHBOR TO SOUTH (6 O'CLOCK) C 1000 - BOUNDARY NEIGHBOR TO WEST (9 O'CLOCK) C THUS J=9 (1001 IN BINARY) IMPLIES THAT THERE ARE BOUNDARY C NEIGHBORS TO THE NORTH AND WEST C EXAMPLES ( X = BNDRY PT., 0 = GRID PT ) C C X X X 0 X C C X 0 0 X X C J=9 J=3 J=14 C C ***** NOTE THAT GTYPE IS INITIALLY SET NEG. AND THEN MADE C POSITIVE WHEN THE INTERIOR IS FILLED C C = 999 C MEANS GRID POINT IS INTERIOR TO THE REGION AND NOT CLOSE C TO THE BOUNDARY, NSIDE HAS BEEN SET TO 999 SO THAT THE PACKING C C = INTEGER LESS THAN 1000 ( 1000 IS IPACKB VALUE) C GRID POINT IS ALSO BOUNDARY PT., GTYPE IS ITS INDEX C C = 0 C GRID POINT IS EXTERIOR FAR FROM THE BOUNDARY C C = NEGATIVE INTEGER C GRID POINT IS EXTERIOR NEXT TO THE BOUNDARY, ITS LOCATION C RELATIVE TO THE BOUNDARY IS ENCODED AS FOR INTERIOR POINTS C NEAR THE BOUNDARY C C 2** BOUNDARY SPECIFICATION ********** C C NBNDPT = NUMBER OF BOUNDARY POINTS ACTUALLY FOUND C XBOUND(I),YBOUND(I) = COORDINATES OF I-TH BOUNDARY POINT C BPARAM(I) = PARAMETER VALUE OF I-TH BOUNDARY POINT C PIECE(I) = INDEX OF BOUNDARY PIECE TO WHICH PT. BELONGS C SMALLEST NUMBER FOR CORNER POINTS C BPTYPE(I) = TYPE OF BOUNDARY POINT C = HORZ,VERT,BOTH,INTE OR JUMP C BNEIGH(I) = POINTER TO THE INTERIOR POINTS FROM THE I-TH C BOUNDARY POINT. SAME SCHEME IS USED TO ENCODE C DIRECTIONS AS FOR THE J PART OF GTYPE ABOVE C BGRID(I) = IX + IPACKB*JY WHEN PT. I IS IN GRID SQUARE IX,JY C C MAXIMUM NUMBER OF BOUNDARY POINTS = NBDIM C THIS IS ESTIMATED BY THE CALLING PROGRAM FOR DIMENSIONING C THE VARIOUS ARRAYS. C C *********** DOMAIN PROCESSOR SUBPROGRAMS *********************** C C REGION - MAIN DRIVER AND USER INTERFACE C HOLE - ALTERNATE DRIVER, INSERTS HOLES IN FIRST DOMAIN C REMOVH - REMOVES HOLE FROM GRID TYPES, UPDATES ALL INFO C DOMAIN - MAIN PROGRAM TO PROCESS A DOMAIN C BWALK - WALK ALONG BOUNDARY TO FIND GRID INTERSECTION C CHKTAN - CHECK FOR BOUNDARY TANGENT TO A GRID LINE C CROSS2 - CHECK FOR DOUBLE CROSSING OF A GRID LINE C DBACK - FOLLOW BOUNDARY FOR A DOUBLE CROSSING OF A GRID LINE C REGULA - MODIFIED REGULA FALSI METHOD C SECANT - SECANT METHOD C CHANGE - MAKE CHANGE OF BOUNDARY PIECE C FILL - LOCATE AND FILL INTERIOR OF THE DOMAIN C EXPAND - EXPAND INTERIOR OF DOMAIN BY 1 POINT C GVALUS - SET TYPES OF ALL GRID POINTS C ISETGT - COMPUTE GTYPE VALUE FOR A GRID POINT, SET IT C INSIDE - UTILITY: CHECKS POINT INSIDE SPECIFIED SUBGRID C LOCATE - UTILITY: LOCATE POINT IN GRID, TYPE IT C NEIGH - COMPUTE POINTERS FROM BOUNDARY TO GRID POINTS C TABLGT - TABLE THE GTYPE VALUES FOR THE GRID C C ******************************************************************** C ************* OUTPUT CONTROL ************************************* C C LEVEL = 0 ***** FATAL ERROR MESSAGES ************************** C BWALK - BOUNDARY GOES OUTSIDE DOMAIN C - UNABLE TO FIND GRID WHERE BOUNDARY GOES C - UNABLE TO FIND GRID INTERSECTION WITH BOUNDARY C CHANGE - BOUNDARY PIECES DO NOT JOIN UP C DOMAIN - BOUNDARY PARAMETER NOT INCREASING C - OVERFLOW OF STORAGE ALLOCATED FOR BOUNDARY POINTS C MUST CHANGE DECLARATIONS IN PROGRAM CALLING REGION C - ABNORMAL EXIT FROM BOUNDARY PROCESSING C PROBABLY CANNOT HAPPEN C - BOUNDARY DOES NOT CLOSE C FILL - FAILURE TO FIND INTERIOR OF DOMAIN C - OVERFLOW IN QUEUE FOR FILLING INTERIOR OF DOMAIN C MUST INCREASE QLIMIT AND RECOMPILE SUBROUTINE FILL C GVALUS - HAVE ILLEGAL POINT TYPE. THINGS ARE REALLY MESSED UP C LOCALE - ASKED TO FIND POINT OUTSIDE DOMAIN C INPUT ERROR OR PROGRAM BUG C REMOVH - HOLE IS TOO CLOSE TO BOUNDARY OF CONTAINING DOMAIN C NEED TO USE FINER GRID C REGION - GRID LINES ARE NOT INCREASING, INPUT ERROR C - LESS THAN 2 GRID LINES IN X OR Y DIRECTION, INPUT ERROR C C LEVEL 1 ***** MINIMAL MESSAGES ********************************* C C REGION - NEW PAGE, DOMAIN PROCESSING STARTS C HOLE - NEW PAGE, HOLE PROCESSING STARTS C DOMAIN - NUMBER OF BOUNDARY/GRID INTERSECTIONS C - FATAL ERROR NOTE( IF PRESENT ) C C LEVEL 2 ***** SUMMARY TRACE ************************************ C C CHANGE - CHANGE OF BOUNDARY PIECE C CHKTAN - POINT REPLACEMENT DUE TO TANGENCY TO GRID LINE C DOMAIN - INITIAL BOUNDARY POINT C - BOUNDARY POINT FOUND C - TABLE OF GRID POINT TYPES( FROM TABLGT ) C - SUMMARY OUTPUT OF BOUNDARY POINT ARRAYS C FILL - WARNING, FIRST TRY TO LOCATE INTERIOR FAILS C - FAILURE INFORMATION ( IF NEEDED ) C REMOVH - TABLE OF GRID POINT TYPES ( FROM TABLGT ) C C LEVEL 3 **** SOME DETAILS OF THE PROGRAM OPERATION ************* C C BWALK - START C - WARNING, POINT FOUND IS NOT IN EXPECTED GRID C - BOUNDARY IS VERTICAL/HORIZONTAL C - FINISH C CHANGE - BOUNDARY PIECE CHANGE DATA C CHKTAN - DATA FOR THE TANGENCY TEST( WHEN IT SUCCEEDS ) C CROSS2 - SIGNS FOR DOUBLE CROSSING CHECK C - SUCCESSFUL RESULTS C DOMAIN - BOUNDARY PIECE CHANGE DATA C - HEADING FOR TYPING OF GRID POINTS C FILL - DATA ON SEARCH FOR INTERIOR POINTS C NEIGH - SUMMARY OF EXECUTION C C LEVEL 4 **** MORE DETAILS ABOUT OPERATION ********************** C C BWALK - EXPANSION OF PARAMETER STEP C DBACK - SUMMARY OF EXECUTION C EXPAND - DETAILS OF IDENTIFICATION OF THE INTERIOR POINTS C GVALUS - DATA USED FOR TYPING GRID POINTS C LOCATE - DATA FOR POINTS LOCATED C NEIGH - DATA FOR NEIGHBORS LOCATED C REGULA - SUMMARY OF EXECUTION C SECANT - START INFORMATION C - FINISH INFORMATION C C LEVEL 5 **** DETAILS OF MANY LOOPS ,ETC. ********************** C C BWALK - DATA ON INITIAL ROUGH GUESSES AT INTERSECTIONS C DBACK - LOOP DETAILS C FILL - DATA FOR LOOP TO FIND FIRST INTERIOR POINT C - DATA ON CLASSIFYING INTERIOR POINTS C GVALUS - DATA ON TYPING GRID POINTS C SECANT - BOUNDS USED ON VARIABLES C - DATA WHEN VARIABLE BOUNDS ARE USED C************************************************* END OF OUTPUT SPECS C EXTERNAL BCOORD INTEGER TYPEC REAL XGRID(NGDIMX), YGRID(NGDIMY), BRANGE(2,NPDIM) DIMENSION GTYPE(NGDIMX,NGDIMY), XBOUND(NBDIM), YBOUND(NBDIM), * BPTYPE(NBDIM), BNEIGH(NBDIM), BGRID(NBDIM), BPARAM(NBDIM), * PIECE(NBDIM) C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C C -------------------------- PROCESS A DOMAIN ------------------------ C C CHECK THAT THE BOUNDARY PARAMETERS ARE ALWAYS INCREASING FATAL = .FALSE. DO 10 IB=N1BND,NBOUND IF (BRANGE(2,IB).GT.BRANGE(1,IB)) GO TO 10 C FATAL ERROR IN BOUNDARY PARAMETERIZATION C IF (LEVEL.GE.0) WRITE (MOUTPT,99999) IB 99999 FORMAT (/5(3H **), 35H FATAL ERROR IN BOUNDARY DEFINITION, /9X, * 19HPARAMETER OF PIECE , I4, 14H IS DECREASING) FATAL = .TRUE. 10 CONTINUE IF (FATAL) GO TO 140 IF (LEVEL.GE.4) WRITE (MOUTPT,99998) CLOCKW, ARC, INHOLE, NGRIDX, * NGRIDY, N1BND, N1BDPT, NBOUND, IPACKB, NSIDE, MOUTPT, QLIMIT, * EPSGRD, EPSTAN 99998 FORMAT (/20X, 45HDATA AND CONSTANTS FOR THIS DOMAIN PROCESSING/3X, * 34HSWITCHES CLOCKW, ARC AND INHOLE = , 3L4/3X, 14HNGRIDX,NGRIDY , * 1H=, 2I4, 31H FIRST BOUNDARY POINT, PIECE = , 2I4, 6H WITH , I3, * 16H BOUNDARY PIECES/3X, 37HCONSTANTS IPACKB,NSIDE,MOUTPT,QLIMIT , * 2H= , 4I6/3X, 37HGEOMETRY TOLERANCES EPSGRD, EPSTAN = , 2E15.5) C LOCATE FIRST BOUNDARY POINT IPIECE = N1BND PARAM = BRANGE(1,IPIECE) CALL BCOORD(PARAM, XB, YB, IPIECE) C CALL LOCATE(XB, YB, IX, JY, TYPE, XGRID, YGRID, NGDIMX, NGDIMY) IF (FATAL) GO TO 140 C SET INFORMATION FOR FIRST BOUNDARY POINT NBNDPT = N1BDPT XBOUND(N1BDPT) = XB YBOUND(N1BDPT) = YB PIECE(N1BDPT) = N1BND BPTYPE(N1BDPT) = TYPE BGRID(N1BDPT) = IX + IPACKB*JY BPARAM(N1BDPT) = PARAM IF ((TYPE.EQ.VERT) .OR. (TYPE.EQ.BOTH)) XBOUND(N1BDPT) = XGRID(IX) IF ((TYPE.EQ.HORZ) .OR. (TYPE.EQ.BOTH)) YBOUND(N1BDPT) = YGRID(JY) C DEBUG IF (LEVEL.GE.2) WRITE (MOUTPT,99997) IPIECE, PARAM, XB, YB, TYPE 99997 FORMAT (5X, 8HON PIECE, I4, 7H PARAM , F10.6, 15H, GIVES COORDS , * 2F10.6, 7H, TYPE , A4) C C ----------------------- START LOOP OVER BOUNDARY PIECES ---------- C FOLLOW EACH PIECE OF THE BOUNDARY AND DETERMINE ALL GRID C LINE INTERSECTIONS AND VARIOUS INFORMATION ASSOCIATED WITH C THE INTERSECTION POINTS( = BOUNDARY POINTS) DO 50 IB=N1BND,NBOUND C C CHECK FOR HITTING END OF BOUNDARY IMMEDIATELY IF (PARAM.GT.BRANGE(2,IPIECE)-EPSGRD) GO TO 40 C C ----------- LOOP WHICH WALKS ALONG THE IB-TH BOUNDARY PIECE ---- DO 30 IDUMB=1,NBDIM C LOCATE THE NEXT GRID LINE INTERSECTION BY WALKING C ALONG THE BOUNDARY TO GRID LINES IXNU,JYNU AND PARAM = PNEXT CALL BWALK(XBOUND(NBNDPT), YBOUND(NBNDPT), IX, JY, PNEXT, * XNEXT, YNEXT, IXNU, JYNU, XGRID, YGRID, NGDIMX, NGDIMY, * BRANGE, NPDIM, BCOORD) IF (FATAL) GO TO 140 C C CHECK FOR HAVING REACHED THE END OF THE BOUNDARY C *** EXIT *** THIS IS THE NORMAL EXIT OF LOOP OVER BOUNDARY PIECE IF (IB.EQ.NBOUND .AND. PNEXT.GT.BRANGE(2,IB)-EPSGRD) GO TO 60 C C CHECK FOR GOING PAST THE END OF THIS BOUNDARY PIECES IF (PNEXT.GE.BRANGE(2,IPIECE)-.5E0*EPSGRD) GO TO 40 C C PUT NEW POINT IN BOUNDARY POINT LIST NBNDPT = NBNDPT + 1 IF (NBNDPT.LE.NBDIM) GO TO 20 C C HAVE EXCEEDED BOUNDARY POINT LIMIT C IF (LEVEL.GE.0) WRITE (MOUTPT,99996) NBDIM 99996 FORMAT (/5(3H **), 35H FATAL ERROR IN PROCESSING DOMAIN, , * 16HESTIMATED NUMBER, I4/9X, 29HOF BOUNDARY-GRID INTERSECTION, * 19H POINTS IS TOO LOW,/9X, 30HTHE REMEDY IS TO INCREASE NBDI, * 2HM ) FATAL = .TRUE. GO TO 140 20 CONTINUE XBOUND(NBNDPT) = XNEXT YBOUND(NBNDPT) = YNEXT PIECE(NBNDPT) = IPIECE BPTYPE(NBNDPT) = TYPE BGRID(NBNDPT) = IXNU + IPACKB*JYNU BPARAM(NBNDPT) = PNEXT IX = IXNU JY = JYNU PARAM = PNEXT IF (LEVEL.GE.2) WRITE (MOUTPT,99995) NBNDPT, XBOUND(NBNDPT), * YBOUND(NBNDPT) 99995 FORMAT (2X, 3(2H--), 22H FOUND BOUNDARY POINT , I3, 7H WITH C, * 11HOORDINATES , 2F10.6) IF (LEVEL.GE.3) WRITE (MOUTPT,99994) IPIECE, TYPE, IX, JY 99994 FORMAT (15X, 9HON PIECE , I2, 9H OF TYPE , A4, 9H IN GRID , * 2I3) C C CHECK TO SEE IF THIS IS A POINT OF TANGENCY. C IF SO, REPLACE THE PREVIOUS POINT IN THE BOUNDARY LIST C BY THIS ONE IF IT IS VERY CLOSE BY CALL CHKTAN(XBOUND, YBOUND, BPARAM, BPTYPE, BGRID, NBDIM, * XGRID, YGRID, NGDIMX, NGDIMY, TYPE) 30 CONTINUE C C ---------------- END OF LOOP OVER POINTS ON 1 BOUNDARY PIECE --- C 40 CONTINUE C C HAVE COME TO END OF PIECE, TAKE CARE OF LAST POINT C THIS ALSO INSERTS THE END POINT IN THE BOUNDARY LIST, C LOCATES THE NEXT POINT ON THE BOUNDARY AND PUTS IT IN C THE LIST TOO. C C DO NOT TRY TO CHANGE THE LAST PIECE IF (IB.EQ.NBOUND) GO TO 60 C CALL CHANGE(BRANGE, NPDIM, XGRID, YGRID, NGDIMX, NGDIMY, NBDIM, * XBOUND, YBOUND, PIECE, BPTYPE, BGRID, BPARAM, BCOORD) IF (FATAL) GO TO 140 CALL LOCATE(XBOUND(NBNDPT), YBOUND(NBNDPT), IX, JY, TYPE, * XGRID, YGRID, NGDIMX, NGDIMY) 50 CONTINUE C C ABNORMAL EXIT FROM BOUNDARY PROCESSING C THIS CANNOT HAPPEN UNLESS THINGS ARE REALLY MESSED UP C IF (LEVEL.GE.0) WRITE (MOUTPT,99993) PARAM, IPIECE, NBOUND, * NBNDPT, IX, JY, NGRIDX, NGRIDY, NPDIM, NBDIM 99993 FORMAT (/5(3H **), 37H ABNORMAL EXIT FROM LOOP IN BOUNDARY , * 25HPROCESSING, = FATAL ERROR/9X, 28HHAPPENS ONLY IF BOUNDARY INF, * 21HO IS REALLY MESSED UP/6X, 8H PARAM =, F12.7/6X, 10H IPIECE,NB, * 48HOUND,NBNDPT,IX,JY,NGRIDX,NGRIDY,NPDIM,NBDIM = /6X, 9I5) FATAL = .TRUE. GO TO 140 C 60 CONTINUE C CHECK CLOSING OF THE BOUNDARY IF (ARC) GO TO 70 CALL BCOORD(BRANGE(2,NBOUND), XLAST, YLAST, IPIECE) DIST = ABS(XLAST-XBOUND(N1BDPT)) + ABS(YLAST-YBOUND(N1BDPT)) IF (DIST.LT.3.*EPSGRD) GO TO 70 C IF (LEVEL.GE.0) WRITE (MOUTPT,99992) 99992 FORMAT (/5(3H **), 37H FATAL ERROR, BOUNDARY DOES NOT CLOSE) FATAL = .TRUE. GO TO 140 C C -------------- END LOOP OVER BOUNDARY PIECES --------------------- 70 CONTINUE C C FIX END OF BOUNDARY, DUPLICATE FIRST POINT AT END C OF LIST, CHECK FOR TANGENCY, PREPARE FOR OUTPUT NBPTF = NBNDPT - N1BDPT + 1 NBDF = NBOUND - N1BND + 1 C C COPY FIRST BOUNDARY POINT IN BOUNDARY LIST AS LAST POINT C ONLY IF THERE IS SPACE FOR IT IF (NBNDPT+1.LT.NBDIM) GO TO 80 C WRITE (MOUTPT,99996) NBDIM FATAL = .TRUE. RETURN C C PRINT WARNING IF LAST PIECE HAS ONLY ONE BOUNDARY POINT 80 IF (PIECE(NBNDPT).EQ.NBOUND) GO TO 90 IF (LEVEL.GE.1) WRITE (MOUTPT,99991) 99991 FORMAT (/5(3H **), 40H WARNING, LAST PIECE HAS ONLY 1 BOUNDARY, * 6H POINT) 90 XBOUND(NBNDPT+1) = XBOUND(N1BDPT) YBOUND(NBNDPT+1) = YBOUND(N1BDPT) PIECE(NBNDPT+1) = NBOUND BPTYPE(NBNDPT+1) = BPTYPE(N1BDPT) BGRID(NBNDPT+1) = BGRID(N1BDPT) BPARAM(NBNDPT+1) = BRANGE(2,NBOUND) C C CHECK FOR TANGENCY AS THE BOUNDARY CLOSES C TEMPORARILY INCREASE BOUNDARY POINT COUNT NBNDPT = NBNDPT + 1 CALL CHKTAN(XBOUND, YBOUND, BPARAM, BPTYPE, BGRID, NBDIM, XGRID, * YGRID, NGDIMX, NGDIMY, TYPE) NBNDPT = NBNDPT - 1 NBNDP1 = NBNDPT + 1 C C IF (LEVEL.GE.3) WRITE (MOUTPT,99990) 99990 FORMAT (//20X, 43HSET GTYPE FOR NEIGHBORS OF BOUNDARY POINTS ) C C FIND BOUNDARY NEIGHBORS GRID POINT INFORMATION C THIS SETS VALUES OF GTYPE FOR GRID POINTS NEXT TO THE BOUNDARY C EXCEPT THAT THE SIGN( = INTERIOR/EXTERIOR FLAG ) IS NOT C SET UNTIL FILL IS CALLED LATER. C SET ALL GRID POINT TYPES = 0 DO 110 IX=1,NGRIDX DO 100 JY=1,NGRIDY GTYPE(IX,JY) = 0 100 CONTINUE 110 CONTINUE CALL GVALUS(XGRID, YGRID, NGDIMX, NGDIMY, GTYPE, NBDIM, XBOUND, * YBOUND, BPTYPE, BGRID) C C NOW FILL IN THE TYPES OF GRID POINTS AWAY FROM THE BOUNDARY C IN THE INTERIOR OF THE REGION. SET GTYPE POS. NEXT TO BOUNDARY C IF (.NOT.ARC) CALL FILL(GTYPE, NGDIMX, NGDIMY, NBDIM, BGRID, * BPTYPE, PIECE, XBOUND, YBOUND) IF (FATAL) GO TO 140 C C GO ALONG THE BOUNDARY AND SET BNEIGH VALUES FOR ALL C THE BOUNDARY POINTS. CALL NEIGH(XGRID, YGRID, NGDIMX, NGDIMY, GTYPE, NBDIM, XBOUND, * YBOUND, BNEIGH, BPTYPE, BGRID) C C SUMMARY OUTPUT IF (LEVEL.GE.1) WRITE (MOUTPT,99989) NBPTF, NBDF, NGRIDX, NGRIDY 99989 FORMAT (/5X, 21HBOUNDARY POINTS FOUND, I12/5X, 15HBOUNDARY PIECES, * 6H FOUND, I12/5X, 9HGRID SIZE, I17, 3H BY, I4/5X, 11HEXECUTION S, * 9HUCCESSFUL) IF (LEVEL.LT.2) GO TO 130 WRITE (MOUTPT,99988) 99988 FORMAT (/3X, 19HBOUNDARY PT XBOUND, 8X, 6HYBOUND, 8X, 6HBPARAM, * 33H PIECE BPTYPE BGRID BNEIGH ) DO 120 IB=N1BDPT,NBNDP1 TYPEC = BPTYPE(IB) WRITE (MOUTPT,99987) IB, XBOUND(IB), YBOUND(IB), BPARAM(IB), * PIECE(IB), TYPEC, BGRID(IB), BNEIGH(IB) 120 CONTINUE 99987 FORMAT (I10, 3F14.9, I5, 4X, A4, I8, I6) 130 CONTINUE C PRINT TABLE OF GTYPE VALUES IF (LEVEL.GT.1) CALL TABLGT(MOUTPT, 1, NGRIDX, NGRIDY, GTYPE) RETURN C C ALL FATAL ERRORS JUMP TO 300 140 CONTINUE C IF (LEVEL.GE.0) WRITE (MOUTPT,99986) 99986 FORMAT (/5(3H **), 23H FATAL ERROR IN DOMAIN,/9X, 12HSTOPS PROCES, * 8HSING RUN) IF (LEVEL.GE.2 .AND. N1BDPT.LT.NBNDPT) WRITE (MOUTPT,99985) * NBPTF, NBDF, NGRIDX, NGRIDY 99985 FORMAT (/5X, 21HBOUNDARY POINTS FOUND, I12/5X, 15HBOUNDARY PIECES, * 6H FOUND, I12/5X, 9HGRID SIZE, I17, 3H BY, I4) IF (LEVEL.LT.2 .OR. N1BDPT.GE.NBNDPT) GO TO 160 DO 150 IB=N1BDPT,NBNDPT TYPEC = BPTYPE(IB) WRITE (MOUTPT,99984) IB, XBOUND(IB), YBOUND(IB), BPARAM(IB), * PIECE(IB), TYPEC, BGRID(IB) 99984 FORMAT (I10, 3F14.9, I5, 4X, A4, I8) 150 CONTINUE 160 CONTINUE FATAL = .TRUE. RETURN END SUBROUTINE BWALK(XSTART, YSTART, IX, JY, PNEW, XNEW, YNEW, INEW, BWA 10 * JNEW, XGRID, YGRID, NGDIMX, NGDIMY, BRANGE, NPDIM, BCOORD) C C THIS SUBROUTINE STARTS AT (XSTART,YSTART) AND WALKS ALONG THE C BOUNDARY TO THE NEXT GRID LINE INTERSECTION AT THE POINT C (XNEW,YNEW) WITH PARAMETER PNEW IN GRID SQUARE INEW,JNEW. C THE PARAMETER START IS PARAM IN COMMON C THE INTERSECTION POINT TYPE IS FOUND AND PLACED IN COMMON AS TYPE C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C EXTERNAL BCOORD LOGICAL XCASE, INSIDE, ALLIN REAL XGRID(NGDIMX), YGRID(NGDIMY), BRANGE(2,NPDIM), GESS1(4), * GESS2(4), TSIGN(4), SIGNOW(4), X(4), Y(4), P(4) INTEGER LTRY(4), TYPEC INTEGER OUTCM(4), WBACK, NOTRY, SUCCES, FAIL, OUTSID DATA WBACK, NOTRY, SUCCES, FAIL, OUTSID /4HBACK,4HNOTR,4HSUCC, * 4HFAIL,4HOUTS/ C C ********** VARIABLES IN BWALK ********** C C *** INPUT * XSTART,YSTART,IX,JY C VIA COMMON = PARAM,EPSGRD,EPSTAN,IPIECE C THE POINT WHERE THE WALK STARTS IS (XSTART,YSTART) C IN GRID SQUARE (IX,JY). C C *** OUTPUT * PNEW = PARAM VALVE FOR NEXT GRID LINE INTERSECTION C INEW, JNEW = INDICES FOR NEXT X,Y GRID LINES C XNEW,YNEW = COORDINATES FOR NEXT INTERSECTION C TYPE(VIA COMMON) = TYPE OF INTERSECTION C = HORZ, VERT OR BOTH C FATAL(VIA COMMON)= SWITCH FOR FATAL ERRORS C C *** MAIN VARIABLES *** C C DELP,DELPNU = CHANGES TRIED FOR PARAM C DXGRD,DYGRD = X AND Y GRID SPACING C XTAR1,XTAR2,YTAR1,YTAR2 = TARGETS OF SECANT METHOD C = COORDINATES OF GRID SQUARE C XTEST,YTEST,PTEST = TEST POINT INSIDE TARGET GRID SQUARE C X(I), 1 TO 4 = X COORDINATES OF SECANT RESULTS C Y(I) 1 TO 4 = Y COORDINATES OF SECANT RESULTS C P(I) 1 TO 4 = PARAMETERS FOR SECANT RESULTS C OUTCM(I) = FLAGS FOR SECANT RESULTS C NOTRY - NOT TRIED C WBACK - WALKED BACK TO START C SUCCES- HIT THE TARGET C FAIL - FAILED TO HIT TARGET C OUTSID- HIT TARGET, BUT OUTSIDE GRID C PEND,XEND,YEND = VALUES OF PARAM,X,Y AT END OF BOUNDARY C PSTOP = PARAMETER VALUE KNOWN TO BE BEYOND C THE NEXT INTERSECTION WITH THE GRID C INEXT, JNEXT = INDEXES OF OTHER TWO SIDES OF SQUARE C BOUNDARY IS PASSING THROUGH C C *** OTHER VARIABLES *** C C NHORZ,NVERT = LOOP COUNTERS C LTRY = ORDER FOR TRIES TO HIT TARGETS C TESTX,TESTY = TEMP VALUES OF X,Y C TSIGN,SIGNOW = SIGNS OF ERROR FOR THE 4 TARGETS C KTAR,STEP = COUNTER AND STEP FOR INITIAL SEARCH C DGMAX,DIST = DISTANCES FROM START C ALLIN = .TRUE. IF ALL PTS MET ARE INSIDE GRID C PEND = BRANGE(2,IPIECE) CALL BCOORD(PEND, XEND, YEND, IPIECE) C SET INITIAL OUTCOME SWITCHES, VARIABLES TO INVALID VALUES DO 10 I=1,4 X(I) = -100.E0 Y(I) = -100.E0 P(I) = PEND + 1.E0 OUTCM(I) = NOTRY 10 CONTINUE ALLIN = .TRUE. C C SET INDICES OF GRIDS THAT MUST CONTAIN NEXT POINT ISOUTH = JY IWEST = IX INORTH = MIN0(JY+1,NGRIDY) IEAST = MIN0(IX+1,NGRIDX) IF (TYPE.EQ.VERT .OR. TYPE.EQ.BOTH) IWEST = MAX0(1,IX-1) IF (TYPE.EQ.HORZ .OR. TYPE.EQ.BOTH) ISOUTH = MAX0(1,JY-1) DXGRD = XGRID(IEAST) - XGRID(IWEST) DYGRD = YGRID(INORTH) - YGRID(ISOUTH) C C FIND PARAMETER CHANGE = PTEST TO GIVE POINT IN THIS GRID FACTOR = 10.E0 PTEST = (PEND-PARAM)/FACTOR DO 30 I1=1,2 DO 20 I2=1,8 CALL BCOORD(PARAM+PTEST, XTEST, YTEST, IPIECE) IF (INSIDE(XTEST,YTEST,IWEST,ISOUTH,IEAST,INORTH,XGRID,NGDIMX, * YGRID,NGDIMY)) GO TO 50 C NOT INSIDE REQUIRED GRID AREA YET C CHECK FOR GOING OUTSIDE GRID DOMAIN IF (.NOT.INSIDE(XTEST,YTEST,1,1,NGRIDX,NGRIDY,XGRID,NGDIMX, * YGRID,NGDIMY)) GO TO 40 ALLIN = .FALSE. PSTOP = PARAM + PTEST PTEST = PTEST/FACTOR 20 CONTINUE FACTOR = 2.E0*FACTOR 30 CONTINUE C C ****** FATAL ERROR ***** C DID NOT FIND REQUIRED POINT = INEXPLICABLE SITUATION C PARAMETER CHANGE IS NOW LESS THAN (PEND-PARAM)/1E+12 C IF (LEVEL.GE.0) WRITE (MOUTPT,99999) XSTART, YSTART 99999 FORMAT (/5(3H **), 35H FATAL ERROR, DID NOT FIND BOUNDARY/9X, * 14HPOINT CLOSE TO, 2F15.9/9X, 31HINEXPLICABLE, IT OCCURS IN FIND, * 23HING INITIAL GRID SQUARE) FATAL = .TRUE. GO TO 520 C C ***** FATAL ERROR, BOUNDARY GOES OUTSIDE GRID ***** 40 CONTINUE C IF (LEVEL.GE.0) WRITE (MOUTPT,99998) XSTART, YSTART 99998 FORMAT (/5(3H **), 27H FATAL ERROR, BOUNDARY GOES/9X, 9HOUTSIDE G, * 13HRID DOMAIN AT, 2F12.6) PTEST = PARAM + PTEST IF (LEVEL.GE.2) WRITE (MOUTPT,99997) XTEST, YTEST, PTEST 99997 FORMAT (9X, 9HTO POINT , 2F12.6, 16H WITH PARAMETER , F12.8) FATAL = .TRUE. GO TO 520 C C NOW CHECK TO SEE IF TEST PT. IS TOO CLOSE TO START 50 DELP2 = PTEST FACTOR = 1.6E0 DO 60 I1=1,6 DX = ABS(XSTART-XTEST) DY = ABS(YSTART-YTEST) IF (AMIN1(DX,DY).GT.4.E0*EPSGRD) GO TO 80 C DO NOT LET DELP2 GET OUT OF RANGE IF (DELP2*FACTOR.GT.PEND-PARAM) GO TO 80 DELP2 = DELP2*FACTOR CALL BCOORD(PARAM+DELP2, XTEST, YTEST, IPIECE) IF (.NOT.INSIDE(XTEST,YTEST,IWEST,ISOUTH,IEAST,INORTH,XGRID, * NGDIMX,YGRID,NGDIMY)) GO TO 70 60 CONTINUE C DID NOT MOVE AWAY FROM START IN BOTH COORDINATES C COULD HAVE VERT OR HORZ BOUNDARY WHICH IS DETECTED LATER GO TO 80 C MOVED TEST POINT OUTSIDE GRID AREA, RESTORE PREVIOUS DELP 70 DELP2 = DELP2/FACTOR ALLIN = .FALSE. 80 PTEST = DELP2 C LOCATE THE POINT INSIDE THE DESIRED GRID AREA CALL BCOORD(PARAM+PTEST, XTEST, YTEST, IPIECE) CALL LOCATE(XTEST, YTEST, IX, JY, TYPEC, XGRID, YGRID, NGDIMX, * NGDIMY) C CHECK LOCATION FAILURE - SHOULD NEVER HAPPEN IF (.NOT.FATAL) GO TO 90 C WRITE (MOUTPT,99996) XSTART, YSTART 99996 FORMAT (/5(3H **), 32H UNABLE TO START SEARCH FOR GRID, 7H INTERS, * 10HECTION AT , 2F12.6) FATAL = .TRUE. GO TO 520 90 INEXT = MIN0(IX+1,NGRIDX) JNEXT = MIN0(JY+1,NGRIDY) C C CHECK IF THE CURRENT BOUNDARY PIECE ENDS INSIDE THE GRID SQUARE IF (INSIDE(XEND,YEND,IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID,NGDIMY) * .AND. ALLIN) GO TO 500 C C HAVE THE GRID SQUARE LOCATED, SET TARGETS XTAR1 = XGRID(INEXT) XTAR2 = XGRID(IX) YTAR1 = YGRID(JNEXT) YTAR2 = YGRID(JY) C COMPUTE THE CORRECT GRID WIDTHS DXGRD = ABS(XGRID(IX)-XGRID(INEXT)) DYGRD = ABS(YGRID(JY)-YGRID(JNEXT)) C IF (LEVEL.GE.3) WRITE (MOUTPT,99995) XSTART, YSTART, IX, JY, * PTEST, XTEST, YTEST, PARAM, XTAR1, XTAR2, YTAR1, YTAR2 99995 FORMAT (5X, 50HBOUNDARY WALK START WITH XSTART,YSTART,IX,JY,PTEST, * 15H,XTEST,YTEST = , 2F11.6, 2I4, 3F11.6/9X, 8HPARAM = , F10.6, * 22H, X AND Y TARGETS ARE , 4F10.6) C C STEP ALONG BOUNDARY TO FIND WHERE IT GOES AND TO OBTAIN ROUGH C ESTIMATES OF INTERSECTIONS. ORDER THE TARGETS WITH MOST C LIKELY FIRST AND OBTAIN INITIAL GUESSES FOR EACH DGMAX = SQRT(DXGRD**2+DYGRD**2) C SET INITIAL SIGNS OF DIRECTION TO TARGETS TSIGN(1) = XTAR1 - XSTART TSIGN(2) = XTAR2 - XSTART TSIGN(3) = YTAR1 - YSTART TSIGN(4) = YTAR2 - YSTART PSTOP = PEND C C **************** RETURN HERE FOR POSSIBLE ADDITIONAL ATTEMPTS C **************** TO BRACKET THE TARGET COORDINATES C LOCAL RETURN FROM INSIDE DO-LOOP 22 C GLOBAL RETURN FROM EXIT OF DO-LOOP 200 100 KTAR = 0 STEP = (PSTOP-PARAM)/5.E0 C SET TARGET ORDERS TO NULL DO 110 I=1,4 LTRY(I) = 0 110 CONTINUE C STEP ALONG THE BOUNDARY DO 160 I=1,5 STEPI = FLOAT(I)*STEP CALL BCOORD(PARAM+STEPI, TESTX, TESTY, IPIECE) DIST = SQRT((TESTX-XSTART)**2+(TESTY-YSTART)**2) IF (DIST.LE.DGMAX) GO TO 120 C HAVE GONE TOO FAR, REDUCE PSTOP PSTOP = PARAM + STEPI C IF THIS IS 1ST OR 2ND POINT, REPEAT WHOLE PROCESS IF (I.LE.2) GO TO 100 C FIND SIGNS TO TARGETS FROM THIS BOUNDARY POINT 120 SIGNOW(1) = XTAR1 - TESTX SIGNOW(2) = XTAR2 - TESTX SIGNOW(3) = YTAR1 - TESTY SIGNOW(4) = YTAR2 - TESTY C C CHECK TO SEE IF ANY SIGN CHANGES HAVE OCCURRED C IF SO, SET ORDER AND SAVE GUESSES DO 150 L=1,4 C SKIP THOSE TARGETS ALREADY IN ORDER IF (KTAR.EQ.0) GO TO 140 DO 130 K=1,KTAR IF (LTRY(K).EQ.L) GO TO 150 130 CONTINUE 140 IF (SIGNOW(L)*TSIGN(L).GE.0.0) GO TO 150 KTAR = KTAR + 1 LTRY(KTAR) = L GESS1(KTAR) = STEPI - STEP GESS2(KTAR) = STEPI 150 CONTINUE 160 CONTINUE C C SET GUESSES AND ORDER FOR TARGETS NOT CROSSED DO 190 L=1,4 C SKIP THOSE TARGETS ALREADY IN ORDER IF (KTAR.EQ.0) GO TO 180 DO 170 K=1,KTAR IF (LTRY(K).EQ.L) GO TO 190 170 CONTINUE 180 KTAR = KTAR + 1 LTRY(KTAR) = L GESS1(KTAR) = 0. GESS2(KTAR) = STEP 190 CONTINUE IF (LEVEL.GE.5) WRITE (MOUTPT,99994) (I,LTRY(I),GESS1(I),GESS2(I), * I=1,4) 99994 FORMAT (9X, 23HTARGET ORDERS,GUESSES =, 4(I3, I2, 2F9.5)) C ----------- CASE STATEMENT TO SEARCH FOR 4 TARGETS ------------ DO 460 I=1,4 KTAR = LTRY(I) GO TO (200, 200, 330, 330), KTAR C C X TARGETS - CHECK FOR VERTICAL BOUNDARY BEFORE SEARCH 200 CONTINUE XCASE = .TRUE. C TEST FOR NEARLY VERTICAL BOUNDARY C C LOOP TO INCREASE DELP AND SEE IF BOUNDARY STAYS C ALONG THE GRID LINE UNTIL DXGRD STEP PASSED DELP = (PSTOP-PARAM)/8.E0 DO 210 NVERT=1,8 CALL BCOORD(PARAM+DELP*FLOAT(NVERT), TESTX, TESTY, IPIECE) DX = ABS(TESTX-XSTART) IF (DX.GT.10.E0*DXGRD*EPSGRD) GO TO 230 IF (DX.GE.DXGRD) GO TO 220 210 CONTINUE C C BOUNDARY IS VERTICAL, SKIP THE X SEARCHES 220 IF (LEVEL.GE.3) WRITE (MOUTPT,99993) XSTART, YSTART 99993 FORMAT (9X, 24HBOUNDARY IS VERTICAL AT , 2F10.6) OUTCM(1) = VERT OUTCM(2) = VERT GO TO 460 C SELECT PROPER X TARGET 230 IF (KTAR.EQ.2) GO TO 270 IF (INEXT.LE.NGRIDX) GO TO 240 C WALKING OUTSIDE GRID DOMAIN IF (LEVEL.GE.1) WRITE (MOUTPT,99992) XSTART, YSTART 99992 FORMAT (5(3H **), 35HPOSSIBLE FATAL ERROR, BOUNDARY GOES, * 24H OUTSIDE GRID DOMAIN AT , 2F11.6) OUTCM(1) = OUTSID OUTCM(2) = OUTSID GO TO 460 240 CONTINUE C CASE 1 -------------------- TRY THE FIRST X-TARGET C SKIP THIS TARGET IF = XSTART, LOOK FOR DOUBLE BACK IF (ABS(XSTART-XTAR1).LE.EPSGRD) GO TO 320 C CHECK FOR TRYING TO WALK OUTSIDE IN X-DIRECTION IF (OUTCM(1).EQ.OUTSID) GO TO 460 C USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS CALL REGULA(XTAR1, XCASE, GESS1(I), GESS2(I), PSTOP, EPSTAN, * BCOORD) C NOW USE THE SECANT METHOD TO HIT TARGET ACCURATELY CALL SECANT(PARAM+GESS1(I), PARAM+GESS2(I), XTAR1, XCASE, P(1), * X(1), Y(1), OUTCM(1), PSTOP, DXGRD, BCOORD) IF (OUTCM(1).EQ.FAIL) GO TO 460 C CHECK FOR HAVING WALKED BACKWARD IF (ABS(XSTART-X(1)).LE.EPSGRD) OUTCM(1) = WBACK C TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE IF (INSIDE(X(1),Y(1),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID, * NGDIMY)) GO TO 250 C POINT FOUND IS OUTSIDE GRID SQUARE, CHECK C FOR HAVING CROSSED THE TARGET TWICE. OUTCM(1) = OUTSID CALL CROSS2(XSTART, XTAR1, XCASE, P(1), X(1), Y(1), OUTCM(1), * BCOORD) C CHECK TO SEE IF THE POINT IS STILL OUTSIDE IF (.NOT.INSIDE(X(1),Y(1),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID, * NGDIMY)) OUTCM(1) = OUTSID 250 IF (OUTCM(1).NE.SUCCES) GO TO 460 C C FOUND CROSSING AT THE NEXT X-GRID LINE 260 PNEW = P(1) XNEW = XTAR1 YNEW = Y(1) JNEW = JY INEW = INEXT C FINISH UP AS IN THE SECOND GRID TARGET GO TO 300 C CASE 2 ---------------- TRY THE SECOND X-TARGET 270 CONTINUE C SKIP THIS TARGET IF = XSTART, LOOK FOR DOUBLE BACK IF (ABS(XSTART-XTAR2).LE.EPSGRD) GO TO 320 C USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS CALL REGULA(XTAR2, XCASE, GESS1(I), GESS2(I), PSTOP, EPSTAN, * BCOORD) C NOW USE THE SECANT METHOD TO HIT TARGET ACCURATELY CALL SECANT(PARAM+GESS1(I), PARAM+GESS2(I), XTAR2, XCASE, P(2), * X(2), Y(2), OUTCM(2), PSTOP, DXGRD, BCOORD) IF (OUTCM(2).EQ.FAIL) GO TO 460 C CHECK FOR HAVING WALKED BACKWARD IF (ABS(XSTART-X(2)).LE.EPSGRD) OUTCM(2) = WBACK C TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE IF (INSIDE(X(2),Y(2),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID, * NGDIMY)) GO TO 280 C POINT FOUND IS OUTSIDE GRID SQUARE, CHECK C FOR HAVING CROSSED THE TARGET TWICE. OUTCM(2) = OUTSID CALL CROSS2(XSTART, XTAR2, XCASE, P(2), X(2), Y(2), OUTCM(2), * BCOORD) C CHECK TO SEE IF THE POINT IS STILL OUTSIDE IF (.NOT.INSIDE(X(2),Y(2),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID, * NGDIMY)) OUTCM(2) = OUTSID 280 IF (OUTCM(2).NE.SUCCES) GO TO 460 C C FOUND CROSSING POINT AT THIS X-GRID LINE 290 CONTINUE JNEW = JY INEW = IX PNEW = P(2) XNEW = XTAR2 YNEW = Y(2) 300 TYPE = VERT C CHECK POSSIBILITY OF TYPE BOTH IF (ABS(YNEW-YTAR1).GT.EPSGRD) GO TO 310 TYPE = BOTH JNEW = JNEXT YNEW = YTAR1 GO TO 510 310 IF (ABS(YNEW-YTAR2).GT.EPSGRD) GO TO 510 TYPE = BOTH JNEW = JY YNEW = YTAR2 GO TO 510 C C SPECIAL CALCULATION FOR DOUBLING BACK ON TARGET 320 CONTINUE CALL DBACK(XSTART, XTEST, PTEST, XCASE, DXGRD, OUTCM(KTAR), * P(KTAR), X(KTAR), Y(KTAR), BCOORD, PSTOP) IF (OUTCM(KTAR).EQ.FAIL) GO TO 460 C CHECK FOR HAVING WALKED BACKWARD TO STARTING POINT IF (ABS(YSTART-Y(KTAR)).LE.EPSGRD) OUTCM(KTAR) = WBACK C TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE IF (.NOT.INSIDE(X(KTAR),Y(KTAR),IX,JY,INEXT,JNEXT,XGRID,NGDIMX, * YGRID,NGDIMY)) OUTCM(KTAR) = OUTSID IF (OUTCM(KTAR).EQ.SUCCES .AND. KTAR.EQ.1) GO TO 260 IF (OUTCM(KTAR).EQ.SUCCES .AND. KTAR.EQ.2) GO TO 290 GO TO 460 C C Y TARGETS - CHECK FOR HORIZONTAL BOUNDARY FIRST 330 CONTINUE XCASE = .FALSE. C C ---------- CHECK FOR NEARLY HORIZONTAL BOUNDARY C LOOP TO INCREASE DELP AND SEE IF BOUNDARY STAYS C ALONG THE GRID LINE UNTIL DYGRD IS PASSED DELP = (PSTOP-PARAM)/8. DO 340 NHORZ=1,8 CALL BCOORD(PARAM+FLOAT(NHORZ)*DELP, TESTX, TESTY, IPIECE) DY = ABS(TESTY-YSTART) IF (DY.GT.10.E0*DYGRD*EPSGRD) GO TO 360 IF (DY.GT.DYGRD) GO TO 350 340 CONTINUE C C BOUNDARY HORIZONTAL 350 IF (LEVEL.GE.3) WRITE (MOUTPT,99991) XSTART, YSTART 99991 FORMAT (9X, 26HBOUNDARY IS HORIZONTAL AT , 2F10.6) OUTCM(3) = HORZ OUTCM(4) = HORZ GO TO 460 C SELECT THE PROPER Y TARGET 360 IF (KTAR.EQ.4) GO TO 400 IF (JNEXT.LE.NGRIDY) GO TO 370 C WALKING OUTSIDE GRID DOMAIN OUTCM(3) = OUTSID OUTCM(4) = OUTSID IF (LEVEL.GE.1) WRITE (MOUTPT,99992) XSTART, YSTART GO TO 460 370 CONTINUE C CASE 3 ------------------ TRY THE FIRST Y-TARGET C SKIP THIS TARGET IF = YSTART, LOOK FOR DOUBLE BACK IF (ABS(YSTART-YTAR1).LE.EPSGRD) GO TO 450 C CHECK FOR TRYING TO WALK OUTSIDE IN Y-DIRECTION IF (OUTCM(3).EQ.OUTSID) GO TO 460 C USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS CALL REGULA(YTAR1, XCASE, GESS1(I), GESS2(I), PSTOP, EPSTAN, * BCOORD) C NOW USE THE SECANT METHOD TO HIT TARGET ACCURATELY CALL SECANT(PARAM+GESS1(I), PARAM+GESS2(I), YTAR1, XCASE, P(3), * X(3), Y(3), OUTCM(3), PSTOP, DYGRD, BCOORD) IF (OUTCM(3).EQ.FAIL) GO TO 460 C CHECK FOR HAVING WALKED BACKWARD IF (ABS(YSTART-Y(3)).LE.EPSGRD) OUTCM(3) = WBACK C TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE IF (INSIDE(X(3),Y(3),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID, * NGDIMY)) GO TO 380 C POINT FOUND IS OUTSIDE GRID SQUARE, CHECK C FOR HAVING CROSSED THE TARGET TWICE. OUTCM(3) = OUTSID CALL CROSS2(YSTART, YTAR1, XCASE, P(3), X(3), Y(3), OUTCM(3), * BCOORD) C CHECK TO SEE IF THE POINT IS STILL OUTSIDE IF (.NOT.INSIDE(X(3),Y(3),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID, * NGDIMY)) OUTCM(3) = OUTSID 380 IF (OUTCM(3).NE.SUCCES) GO TO 460 C C FOUND CROSSING POINT AT THE NEXT Y-GRID LINE 390 CONTINUE JNEW = JNEXT INEW = IX PNEW = P(3) XNEW = X(3) YNEW = YTAR1 C FINISH UP AS AT THE FOURTH TARGET GO TO 430 C CASE 4 ------------ TRY THE SECOND Y-TARGET 400 CONTINUE C SKIP THIS TARGET IF = YSTART, LOOK FOR DOUBLE BACK IF (ABS(YSTART-YTAR2).LE.EPSGRD) GO TO 450 C USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS CALL REGULA(YTAR2, XCASE, GESS1(I), GESS2(I), PSTOP, EPSTAN, * BCOORD) C NOW USE THE SECANT METHOD TO HIT TARGET ACCURATELY CALL SECANT(PARAM+GESS1(I), PARAM+GESS2(I), YTAR2, XCASE, P(4), * X(4), Y(4), OUTCM(4), PSTOP, DYGRD, BCOORD) IF (OUTCM(4).EQ.FAIL) GO TO 460 C CHECK FOR HAVING WALKED BACKWARD IF (ABS(YSTART-Y(4)).LE.EPSGRD) OUTCM(4) = WBACK C TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE IF (INSIDE(X(4),Y(4),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID, * NGDIMY)) GO TO 410 C POINT FOUND IS OUTSIDE GRID SQUARE, CHECK C FOR HAVING CROSSED THE TARGET TWICE. OUTCM(4) = OUTSID CALL CROSS2(YSTART, YTAR2, XCASE, P(4), X(4), Y(4), OUTCM(4), * BCOORD) C CHECK TO SEE IF THE POINT IS STILL OUTSIDE IF (.NOT.INSIDE(X(4),Y(4),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID, * NGDIMY)) OUTCM(4) = OUTSID 410 IF (OUTCM(4).NE.SUCCES) GO TO 460 C C FOUND CROSSING POINT AT THIS Y-GRID LINE 420 CONTINUE PNEW = P(4) XNEW = X(4) YNEW = YTAR2 INEW = IX JNEW = JY 430 TYPE = HORZ C CHECK POSSIBILITY OF TYPE BOTH C WE GET HERE ONLY IN UNUSUAL CASES IF (ABS(XNEW-XTAR1).GT.EPSGRD) GO TO 440 TYPE = BOTH INEW = INEXT XNEW = XTAR1 GO TO 510 440 IF (ABS(XNEW-XTAR2).GT.EPSGRD) GO TO 510 TYPE = BOTH INEW = IX XNEW = XTAR2 GO TO 510 C C SPECIAL CALCULATION FOR DOUBLING BACK ON TARGET 450 CONTINUE CALL DBACK(YSTART, YTEST, PTEST, XCASE, DYGRD, OUTCM(KTAR), * P(KTAR), X(KTAR), Y(KTAR), BCOORD, PSTOP) IF (OUTCM(KTAR).EQ.FAIL) GO TO 460 C CHECK FOR HAVING WALKED BACKWARD TO STARTING POINT IF (ABS(XSTART-X(KTAR)).LE.EPSGRD) OUTCM(KTAR) = WBACK C TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE IF (.NOT.INSIDE(X(KTAR),Y(KTAR),IX,JY,INEXT,JNEXT,XGRID,NGDIMX, * YGRID,NGDIMY)) OUTCM(KTAR) = OUTSID IF (OUTCM(KTAR).EQ.SUCCES .AND. KTAR.EQ.3) GO TO 390 IF (OUTCM(KTAR).EQ.SUCCES .AND. KTAR.EQ.4) GO TO 420 C 460 CONTINUE C C HAVE CHECKED ALL POSSIBLE GRID INTERSECTIONS - NOTHING FOUND C CHECK THE OUTSIDE POINTS FOR ACCEPTABILITY C OR FOR REDUCING PSTOP AND MAKING THE WHOLE SEARCH OVER C AGAIN - GOING BACK UP TO 12 POUT = PSTOP DO 470 I1=1,4 IF (OUTCM(I1).NE.OUTSID) GO TO 470 POUT = AMIN1(P(I1),POUT) IF (.NOT.INSIDE(X(I1),Y(I1),IWEST,ISOUTH,IEAST,INORTH,XGRID, * NGDIMX,YGRID,NGDIMY)) GO TO 470 C HAVE ACCEPTABLE POINT, TAKE IT IF (LEVEL.GE.3) WRITE (MOUTPT,99990) X(I1), Y(I1) 99990 FORMAT (1X, 8(1H+), 10HTHE POINT , 2F10.6, 17H TAKEN EVEN THOUG, * 9HH OUTSIDE, 14H EXPECTED GRID) C LOCATE NEW POINT,SET COORDINATES AND EXIT BWALK CALL LOCATE(X(I1), Y(I1), INEW, JNEW, TYPE, XGRID, YGRID, * NGDIMX, NGDIMY) C CHECK LOCATION FAILURE - SHOULD NEVER HAPPEN IF (FATAL) GO TO 470 XNEW = X(I1) YNEW = Y(I1) PNEW = P(I1) GO TO 510 470 CONTINUE C C CHECK TO SEE IF PSTOP HAS BEEN REDUCED. IF SO REDO THE C WHOLE SEARCH - GOING BACK UP TO 12 TO RESTART IF (POUT.GE.PSTOP) GO TO 490 PSTOP = POUT IF (LEVEL.GE.3) WRITE (MOUTPT,99989) PSTOP, OUTCM(1), XTAR1, * X(1), Y(1), P(1), OUTCM(2), XTAR2, X(2), Y(2), P(2), OUTCM(3), * YTAR1, X(3), Y(3), P(3), OUTCM(4), YTAR2, X(4), Y(4), P(4) 99989 FORMAT (9X, 42HALL POINTS FOUND ARE TOO FAR AWAY. RESTART, * 36H SEARCH WITH SMALLER PARAMETER LIMIT, F11.6/3X, 10HTHE RESULT, * 31HS OF 4 TRIES TO FOLLOW BOUNDARY/2(2X, A4, 4H X=, F10.6, * 2F14.8, F12.6)/2(2X, A4, 4H Y=, F10.6, 2F14.8, F12.6)) DO 480 I=1,4 X(I) = -100.E0 Y(I) = -100.E0 P(I) = PEND + 1.E0 OUTCM(I) = NOTRY 480 CONTINUE ALLIN = .TRUE. GO TO 100 C STILL NO POINT FOUND C ************************************* C *********** FATAL ERROR *********** C ************************************* 490 CONTINUE C IF (LEVEL.GE.0) WRITE (MOUTPT,99988) IPIECE, XSTART, YSTART, IX, * JY, OUTCM(1), XTAR1, X(1), Y(1), P(1), OUTCM(2), XTAR2, X(2), * Y(2), P(2), OUTCM(3), YTAR1, X(3), Y(3), P(3), OUTCM(4), YTAR2, * X(4), Y(4), P(4) 99988 FORMAT (/, 5(3H **), 24H FATAL ERROR IN BOUNDARY, 1X, /, 3X, * 29HOR GRID DEFINITION FOR PIECE , I3, 10H AT POINT , 2F10.6, /, * 2X, 9H IN GRID , 2I3, 3X, 32HTHE RESULTS OF 4 TRIES TO FOLLOW, * 9H BOUNDARY, /, 2X, 41HNOTR - TARGET NOT TRIED OUTS - TAR, * 17HGET OUTSIDE GRID,/, 2X, 34HFAIL - TARGET NOT ACHIEVED BAC, * 24HK - FOUND PREVIOUS PT. ,/, 2X, 7HOUTCOME, 4X, 6HTARGET, 12X, * 10HX,Y-COORDS, 10X, 9HPARAMETER, /, 2(2X, A4, 4H X=, F10.6, * 2F14.8, F12.6, /), 2(2X, A4, 4H Y=, F10.6, 2F14.8, F12.6, /)) FATAL = .TRUE. GO TO 520 C C CURRENT PIECE STOPS IN THE GRID SQUARE 500 CONTINUE TYPE = INTER IF (LEVEL.GE.2) WRITE (MOUTPT,99987) IPIECE, XEND, YEND, IX, JY 99987 FORMAT (2X, 3(2H--), 6H PIECE, I3, 14H ENDS AT POINT, 2F10.6, * 12H NEAR GRID , 2I4) XNEW = XEND YNEW = YEND PNEW = PEND + .1E0*EPSGRD C C NORMAL SUBPROGRAM EXIT 510 CONTINUE C PRINT POSSIBLE BEFORE RETURN IF (LEVEL.GE.3) WRITE (MOUTPT,99986) PNEW, XNEW, YNEW, TYPE, * INEW, JNEW, (OUTCM(I),X(I),Y(I),I=1,4) 99986 FORMAT (5X, 32H--- END OF BOUNDARY WALK, PNEW =, F12.7/5X, * 12HGIVES POINT , 2F10.6, 9H OF TYPE , A4, 9H IN GRID , 2I3, * 27H OUTCOMES,PTS FOR 4 CASES =/3X, 4(2X, A4, 2F12.6)) RETURN C EXIT FOR FATAL ERROR DETECTED AT 200 ABOVE 520 FATAL = .TRUE. RETURN END SUBROUTINE CHANGE(BRANGE, NPDIM, XGRID, YGRID, NGDIMX, NGDIMY, CHA 10 * NBDIM, XBOUND, YBOUND, PIECE, BPTYPE, BGRID, BPARAM, BCOORD) C C *** THIS PROGRAM CHANGES BOUNDARY PIECES AND CHECKS FOR CONTINUTIY C IT ALSO PUTS THE END OF THE BOUNDARY PLUS THE NEXT POINT IN THE C BOUNDARY LIST. C C INPUT = BRANGE, (VIA COMMON) IPIECE,PARAM C PLUS DIMENSION VARIABLES = NPDIM,NGDIMX,NGDIMY,NBDIM C OUTPUT= (VIA COMMON) TYPE,PARAM,IPIECE,NBNDPT C XBOUND(NBNDPT) PIECE(NBNDPT) BPTYPE(NBNDPT) C YBOUND(NBNDPT) BGRID(NBNDPT) BPARAM(NBNDPT) C C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C EXTERNAL BCOORD INTEGER TYPEC REAL XGRID(NGDIMX), YGRID(NGDIMY), BRANGE(2,NPDIM) DIMENSION XBOUND(NBDIM), YBOUND(NBDIM), PIECE(NBDIM), * BPTYPE(NBDIM), BGRID(NBDIM), BPARAM(NBDIM) C COORDINATES OF END OF CURRENT PIECE CALL BCOORD(BRANGE(2,IPIECE), XOLD, YOLD, IPIECE) C C INCREMENT COUNTER IPIECE OF BOUNDARY PIECES C FIND COORDINATES (XCOR,YCOR) THAT START THE NEW PIECE IPIECE = IPIECE + 1 PARAM = BRANGE(1,IPIECE) CALL BCOORD(PARAM, XCOR, YCOR, IPIECE) C C CHECK THAT THE NEXT PIECE JOINS UP WITH CURRENT PIECE DIST = SQRT((XOLD-XCOR)**2+(YOLD-YCOR)**2) IF (DIST.LE.4.E0*EPSGRD) GO TO 10 C IF (LEVEL.GE.0) WRITE (MOUTPT,99999) IPIECE, XCOR, YCOR, DIST 99999 FORMAT (/5(3H **), 20H FATAL ERROR, PIECE , I3, 10H STARTS AT, * 2F15.7/9X, 32HPIECES DO NOT JOIN UP, BREAK IS , 1PE15.6) FATAL = .TRUE. RETURN C 10 CONTINUE C C LOCATE AND TYPE CORNER POINT CALL LOCATE(XCOR, YCOR, IXCOR, JYCOR, TYPEC, XGRID, YGRID, * NGDIMX, NGDIMY) C C DEBUG IF (LEVEL.GE.3) WRITE (MOUTPT,99998) IPIECE, NBNDPT, PARAM, XCOR, * YCOR, TYPEC, IXCOR, JYCOR 99998 FORMAT (1X, 18H---CHANGE TO PIECE, I3, 18H AT BOUNDARY POINT, I3, * 14H WITH PARAM = , F10.6, 18H, CORNER POINT IS , 2F10.6, * 10H, OF TYPE , A4/4X, 14HIT IS IN GRID , 2I3) C C CHECK TO SEE IF THE END POINT IS ALREADY IN BOUNDARY LIST C IF SO, GO ON TO THE NEXT POINT IF (ABS(BRANGE(2,IPIECE-1)-BPARAM(NBNDPT)).LE.EPSGRD) GO TO 20 C C HAVE NEW POINT (XCOR,YCOR) TO INSERT IN BOUNDARY LIST NBNDPT = NBNDPT + 1 XBOUND(NBNDPT) = XCOR YBOUND(NBNDPT) = YCOR BPTYPE(NBNDPT) = TYPEC BGRID(NBNDPT) = IXCOR + IPACKB*JYCOR BPARAM(NBNDPT) = BRANGE(2,IPIECE-1) PIECE(NBNDPT) = IPIECE - 1 C CHECK FOR TANGENCY AT CORNER POINT AS CHANGE IS MADE CALL CHKTAN(XBOUND, YBOUND, BPARAM, BPTYPE, BGRID, NBDIM, XGRID, * YGRID, NGDIMX, NGDIMY, TYPEC) C C ******* SUBPROGRAM EXIT AND DEBUG OUTPUT ******* 20 CONTINUE IF (LEVEL.LE.1) RETURN TYPEC = BPTYPE(NBNDPT) WRITE (MOUTPT,99997) IPIECE, XCOR, YCOR, TYPEC 99997 FORMAT (/9X, 25HCHANGE TO BOUNDARY PIECE , I3, 4H AT , 2F10.6, * 13H WITH TYPE = , A4) WRITE (MOUTPT,99996) NBNDPT, XBOUND(NBNDPT), YBOUND(NBNDPT) 99996 FORMAT (2X, 3(2H--), 22H IT IS BOUNDARY POINT , I3, 10H WITH COOR, * 8HDINATES , 2F10.6) IF (LEVEL.EQ.2) RETURN WRITE (MOUTPT,99995) IXCOR, JYCOR, TYPEC, BPARAM(NBNDPT), DIST 99995 FORMAT (5X, 43HCORNER POINT DATA IX,JY,TYPE,PARAM,DISTANCE, 2I4, * 1X, A4, F10.6, 1PE10.2) RETURN END SUBROUTINE CHKTAN(XBOUND, YBOUND, BPARAM, BPTYPE, BGRID, NBDIM, CHK 10 * XGRID, YGRID, NGDIMX, NGDIMY, TYPEPT) C C ***** THIS PROGRAM CHECKS FOR THE BOUNDARY BEING TANGENT TO A GRID C LINE. IF TWO BOUNDARY POINTS ARE VERY CLOSE, THAT IS ABOUT C SQRT(EPSGRD) APART. THEN THEY ARE ASSUMED TO BE THE SAME POINT C AND THE FIRST POINT IS DISCARDED AND REPLACED BY THE LAST ONE C FOUND( HAS INDEX = NBNDPT ). C REAL XBOUND(NBDIM), YBOUND(NBDIM), BPARAM(NBDIM), XGRID(NGDIMX), * YGRID(NGDIMY) DIMENSION BPTYPE(NBDIM), BGRID(NBDIM) INTEGER TYPEPT, TYPEC INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C DELX = ABS(XBOUND(NBNDPT)-XBOUND(NBNDPT-1)) DELY = ABS(YBOUND(NBNDPT)-YBOUND(NBNDPT-1)) C CHECK FOR TANGENCY TO HORIZONTAL GRID LINE C IF POINT IS ON HORIZONTAL GRID LINE IF (.NOT.(TYPEPT.EQ.HORZ .OR. TYPEPT.EQ.BOTH)) GO TO 10 IF (DELX.LE.EPSTAN .AND. DELY.LE.EPSGRD) GO TO 30 C C CHECK FOR TANGENCY TO VERTICAL GRID LINE C IF POINT IS ON VERTICAL GRID LINE 10 IF (.NOT.(TYPEPT.EQ.VERT .OR. TYPEPT.EQ.BOTH)) GO TO 20 IF (DELY.LE.EPSTAN .AND. DELX.LE.EPSGRD) GO TO 30 C C NO TANGENCY, SUBROUTINE EXIT 20 RETURN C C HAVE FOUND TANGENCY, MERGE THE 2 POINTS C ITEMS NOT CHANGED HERE ARE THE SAME FOR BOTH POINTS 30 CONTINUE XBOUND(NBNDPT-1) = XBOUND(NBNDPT) YBOUND(NBNDPT-1) = YBOUND(NBNDPT) BPTYPE(NBNDPT-1) = BPTYPE(NBNDPT) BPARAM(NBNDPT-1) = BPARAM(NBNDPT) BGRID(NBNDPT-1) = BGRID(NBNDPT) NBNDPT = NBNDPT - 1 TYPEC = BPTYPE(NBNDPT) IF (LEVEL.GE.2) WRITE (MOUTPT,99999) NBNDPT, XBOUND(NBNDPT), * YBOUND(NBNDPT), TYPEC 99999 FORMAT (2X, 4(3H+++), 40H FOUND BOUNDARY TANGENT TO GRID AT POINT, * I4/9X, 2H= , 2F10.6, 9H OF TYPE , A4, 23H REPLACES PREVIOUS POIN, * 1HT) IF (LEVEL.GE.3) WRITE (MOUTPT,99998) DELX, DELY 99998 FORMAT (8X, 27H TEST VARIABLES DELX,DELY =, 4E12.3) RETURN END SUBROUTINE CROSS2(START, TARGET, XCASE, PEND, XEND, YEND, OUTCOM, CRO 10 * BCOORD) C C **** THIS ROUTINE IS TO CHECK FOR THE POSSIBILITY THAT THE BOUNDARY C CROSSED THE TARGET GRID LINE MORE THAT ONCE AND THAT THE C FIRST CROSSING WAS NOT LOCATED. THE IDEA IS TO CHECK TO SEE C IF THE BOUNDARY IS MOVING TOWARD OR AWAY FROM THE START. IF C A DOUBLE( OR MORE ) CROSSING IS DETECTED THEN REGULA ( THE C MODIFIED REGULA FALSI METHOD) IS USED TO LOCATE A CROSSING C CLOSER TO THE START. NOTE THIS LOGIC WILL NOT NOT DETECT A C TRIPLE CROSSING. C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP INTEGER OUTCOM, SUCCES EXTERNAL BCOORD LOGICAL XCASE DATA SUCCES /4HSUCC/ C C FIRST CHECK FOR A DOUBLE CROSSING PCHECK = PEND - EPSTAN CALL BCOORD(PCHECK, XCHECK, YCHECK, IPIECE) CHECK = YCHECK IF (XCASE) CHECK = XCHECK SIGN = TARGET - START CSIGN = TARGET - CHECK IF (LEVEL.GE.3) WRITE (MOUTPT,99999) PEND, SIGN, CSIGN 99999 FORMAT (9X, 6H***** , 24HDOUBLE CROSSING CHECK AT, F12.6, * 10H HAS SIGNS, 2E13.3) IF (SIGN*CSIGN.GT.0.) RETURN C FOUND A SIGN CHANGE = THERE IS CLOSER CROSSING C USE REGULA TO FIND IT ACCURATELY PLEFT = 0. PRIGHT = PCHECK - PARAM CALL REGULA(TARGET, XCASE, PLEFT, PRIGHT, PCHECK-PARAM, EPSGRD, * BCOORD) C REGULA SHOULD ALWAYS CONVERGE SINCE THE TARGET IS BRACKETED C CHECK CONVERGENCE CALL BCOORD(PARAM+PRIGHT, X, Y, IPIECE) V = Y IF (XCASE) V = X IF (ABS(TARGET-V).GT.EPSGRD) RETURN C C HAVE CONVERGENCE, SET VALUES AND OUTCOM OUTCOM = SUCCES IF (LEVEL.GE.3) WRITE (MOUTPT,99998) XEND, YEND, X, Y 99998 FORMAT (9X, 6H***** , 40HDOUBLE CROSSING CHECK FOUND CLOSER POINT/ * 15X, 17HORIG OUTSIDE PT =, 2F12.6, 9H NEW PT =, 2F12.6) PEND = PRIGHT XEND = X YEND = Y RETURN END SUBROUTINE DBACK(TARGET, VTEST, DELPST, XCASE, DVGRD, OUTCOM, P, DBA 10 * X, Y, BCOORD, PSTOP) C C ***** THIS PROGRAM ATTEMPTS TO FIND A CROSSING DUE TO THE C BOUNDARY DOUBLING BACK OVER A GRID LINE WITHIN THE GRID C SQUARE. THE BOUNDARY IS FOLLOWED FROM (XTEST,YTEST) UNTIL C IT CROSSES THE TARGET AGAIN. THIS GIVES PARAMETERS P1,P2 C WHICH ARE USED TO START SECANT SEARCH FOR THE RESULTS C P,X,Y = SOLUTION. C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C EXTERNAL BCOORD INTEGER OUTCOM, FAIL LOGICAL XCASE DATA FAIL /4HFAIL/ C ALWAYS LIMIT PARAMETER TO CURRENT PIECE DELPMX = PSTOP - PARAM FACTOR = 1.5E0 C DEVST = SIZE OF POINT-TARGET AT START C DEVNU = SIZE OF POINT-TARGET AT NEXT STEP DEVST = VTEST - TARGET DELP = DELPST KEXPND = 0 IF (LEVEL.GE.5) WRITE (MOUTPT,99999) KEXPND, DELP, DEVST, VTEST C C EXPAND PARAMETER STEP UNTIL BOUNDARY C CROSSES THE TARGET GRID LINE DO 10 KEXPND=1,20 C STEP AWAY FROM PREVIOUS POINT TO CROSS THE TARGET DELPNU = AMIN1(FACTOR*DELP,DELPMX) CALL BCOORD(PARAM+DELPNU, XNU, YNU, IPIECE) VNU = YNU IF (XCASE) VNU = XNU DEVNU = VNU - TARGET C TEST FOR FINDING A SIGN CHANGE IN POINT-TARGET DIFFERENCE C THIS MEANS THAT WE HAVE CROSSED THE TARGET GRID LINE C SO WE START TO USE SECANT METHOD IF (DEVST*DEVNU.LT.0.E0) GO TO 30 DELP = DELPNU C INCREASE FACTOR OCCAISIONALLY IF (MOD(KEXPND,4).EQ.0) FACTOR = FACTOR*1.2E0 C CHECK FOR GETTING TOO FAR AWAY = FAILURE IF (ABS(DEVNU).GE.DVGRD .OR. DELPNU.GE.DELPMX) GO TO 20 IF (LEVEL.GE.5) WRITE (MOUTPT,99999) KEXPND, DELPNU, DEVNU, * XNU, YNU 10 CONTINUE 99999 FORMAT (9X, 10HDBACK LOOP, I3, 25H GIVES DELPNU,DEVNU,X,Y =, * 4F11.6) C C FAILED TO FIND A DOUBLE CROSSING = FAILURE EXIT 20 OUTCOM = FAIL IF (LEVEL.GE.4) WRITE (MOUTPT,99998) TARGET, VTEST, XNU, YNU 99998 FORMAT (9X, 32HDBACK ROUTINE FAILS FOR TARGET , F10.6, 7H, START, * 4H AT , F11.7, 9H, END AT , 2F11.7) RETURN C C FOUND CROSSING, SET PARAMETER GUESSES P1 AND P2 C CALL SECANT TO FINISH FINDING CROSSING 30 CONTINUE P1 = PARAM + DELPNU P2 = PARAM + DELP IF (LEVEL.GE.5) WRITE (MOUTPT,99997) P1, P2 99997 FORMAT (9X, 33HDBACK SUCCEEDED TO BRACKET TARGET, 12H WITH PARAME, * 6HTER = , 2F10.6) C USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS CALL REGULA(TARGET, XCASE, DELPNU, DELP, PSTOP, EPSTAN, BCOORD) C USE SECANT TO GET THE FINAL RESULTS CALL SECANT(PARAM+DELPNU, PARAM+DELP, TARGET, XCASE, P, X, Y, * OUTCOM, PSTOP, DVGRD, BCOORD) IF (OUTCOM.EQ.FAIL .AND. LEVEL.GE.4) WRITE (MOUTPT,99998) TARGET, * VTEST, XNU, YNU RETURN END SUBROUTINE EXPAND(ICENT, JCENT, NUPTS, INEW, JNEW, GTYPE, NGDIMX, EXP 10 * NGDIMY) C C ***** THIS PROGRAM IS GIVEN THE POINT ICENT,JCENT AS AN INTERIOR C POINT. IT THEN DETERMINES IF ITS 4 NEIGHBORS ARE INTERIOR C POINTS AND, IF SO, MARKS THEM BY C GTYPE = NSIDE FOR PTS AWAY FROM THE BOUNDARY C GTYPE = -GTYPE FOR PTS NEXT TO THE BOUNDARY C C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C DIMENSION GTYPE(NGDIMX,NGDIMY) C INTEGER INEW(4), JNEW(4), PTYPE, IXN(4), JYN(4), CENTER LOGICAL PTNAY, BCROSS DATA IXN(1), IXN(2), IXN(3), IXN(4) /0,-1,0,1/ DATA JYN(1), JYN(2), JYN(3), JYN(4) /-1,0,1,0/ C C ARITHMETIC STATEMENT FUNCTIONS C PTNAY IS TRUE FOR BOUNDARY POINT LOCATED C ABOVE J=1 C RIGHT J=2 C BELOW J=3 C LEFT J=4 PTNAY(K,J) = MOD(K/(IPACKB*2**(MOD(J-1,4))),2).EQ.1 C C BCROSS IS TRUE IF THE LINE FROM THE J-TH POINT TO THE C CENTER POINT CROSSES THE BOUNDARY BCROSS(J) = PTNAY(PTYPE,J) .AND. PTNAY(CENTER,J+2) C CENTER = GTYPE(ICENT,JCENT) NUPTS = 0 DO 10 JNAY=1,4 IX = ICENT + IXN(JNAY) JY = JCENT + JYN(JNAY) PTYPE = -GTYPE(IX,JY) C PTYPE IS NEGATIVE FOR BOUNDARY POINTS IF (PTYPE.LT.0) GO TO 10 IF (BCROSS(JNAY)) GO TO 10 C HAVE AN INTERIOR POINT C CHECK FOR LEGALITY OF INDEXES IF (IX.LT.1 .OR. IX.GT.NGRIDX) GO TO 10 IF (JY.LT.1 .OR. JY.GT.NGRIDY) GO TO 10 NUPTS = NUPTS + 1 INEW(NUPTS) = IX JNEW(NUPTS) = JY IF (PTYPE.EQ.0) PTYPE = NSIDE GTYPE(IX,JY) = PTYPE 10 CONTINUE IF (LEVEL.GE.4) WRITE (MOUTPT,99999) ICENT, JCENT, PTYPE 99999 FORMAT (5X, 19HEXPAND INTERIOR AT , 2I3, 13H FINAL TYPE =, I7) IF (LEVEL.GE.5 .AND. NUPTS.GT.0) WRITE (MOUTPT,99998) * (I,INEW(I),JNEW(I),I=1,NUPTS) 99998 FORMAT (9X, 14HNEW POINTS = , 4(I7, 3H = , 2I3)) RETURN END SUBROUTINE FILL(GTYPE, NGDIMX, NGDIMY, NBDIM, BGRID, BPTYPE, FIL 10 * PIECE, XBOUND, YBOUND) C C ***** THIS PROGRAM LOCATES AN INTERIOR POINT NEAR THE START OF THE C BOUNDARY AND THEN EXPANDS THE INTERIOR, MARKING ALL INTERIOR C POINTS BY CHANGING GTYPE TO PLUS OR SETTING GTYPE = INTER. C THE METHOD USES A QUEUE OF POINTS ON THE FRINGE OF THE CURRENTLY C KNOWN INTERIOR AND THIS QUEUE IS PROCESSED VIA EXPAND WHICH C LOCATES NEW INTERIOR POINTS ADJACENT TO ANY GIVEN ONE. C C INPUT. XBOUND,YBOUND = BOUNDARY POINT COORDINATES C BGRID = - - GRID LOCATIONS C BPTYPE = - - TYPES C PIECE = - - PIECE C GTYPE = GRID POINT MARKERS C OUTPUT. GTYPE MODIFIED C DIMENSION GTYPE(NGDIMX,NGDIMY) DIMENSION BGRID(NBDIM), BPTYPE(NBDIM), XBOUND(NBDIM), * YBOUND(NBDIM), PIECE(NBDIM) C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C INTEGER INEW(4), JNEW(4) INTEGER TYPE1, TYPE2, TYPE3, CORNER LOGICAL FIRST DATA CORNER /4HCORN/ C C SET ORIENT = +- 1 ( A SIGN ) DEPENDING ON ORIENTATION ORIENT = -1. IF (CLOCKW) ORIENT = 1. C CHANGE ORIENTATION IF CALLED FROM HOLE IF (INHOLE) ORIENT = -ORIENT C C FIND AN INITIAL INTERIOR POINT FROM THE BPTYPE ARRAY IB = N1BDPT FIRST = .TRUE. C C RETURN HERE FOR RESTART OF SEARCH IF THINGS LOOK C BAD. SKIP TO END OF LOOP IF IB=NBNDPT. 10 CONTINUE ISTART = IB + 1 IF (ISTART.LT.NBNDPT) GO TO 30 C IF THIS IS THE SECOND PASS OVER THE BOUNDARY, QUIT IF (.NOT.FIRST) GO TO 60 C HAVE FAILED TO FIND INTERIOR POINT. MAKE ANOTHER PASS C AROUND BOUNDARY NOW ALLOWING CORNER POINTS TO BE USED C NEXT TO THE BOUNDARY POINT BEING TESTED. 20 IF (LEVEL.GE.2) WRITE (MOUTPT,99999) 99999 FORMAT (2X, 7(3H++ ), 36HFIRST ATTEMPT TO FIND INTERIOR POINT, * 7H FAILED) FIRST = .FALSE. ISTART = N1BDPT + 1 30 CONTINUE C DEBUG FOR LOOP RESTART IF (LEVEL.GE.4 .AND. ISTART.GT.N1BDPT+1) WRITE (MOUTPT,99998) * IADDX, TESTX, IADDY, TESTY, IDIFF 99998 FORMAT (4X, 44H-RESTART SEARCH WITH IADDX,TESTX,IADDY,TESTY, * 8H,IDIFF =, 2(I3, E14.3), I4) C C USE TYPES FOR 3 POINTS IN A ROW, INITIALIZE 2 HERE C SET TYPE FOR CORNER POINTS = CORNER TYPE2 = CORNER TYPE3 = BPTYPE(ISTART) IF (PIECE(ISTART).NE.PIECE(ISTART+1)) TYPE3 = CORNER DO 50 IB=ISTART,NBNDPT C UPDATE TYPES OF 3 POINTS ABOUT IB POINT TYPE1 = TYPE2 TYPE2 = TYPE3 TYPE3 = BPTYPE(IB+1) IF (PIECE(IB+1).NE.PIECE(IB+2)) TYPE3 = CORNER C DO NOT LOOK FOR INTERIOR POINT NEXT TO BOUNDARY POINT C WHOSE SUCCESSOR OR PREDECESSOR IS A CORNER C I. E. WHOSE TYPE IS CORNER OR INTER. IF (TYPE2.EQ.CORNER .OR. TYPE2.EQ.INTER) GO TO 40 IF (FIRST .AND. (TYPE3.EQ.CORNER .OR. TYPE1.EQ.CORNER)) GO TO 40 IF (TYPE3.EQ.INTER .OR. TYPE1.EQ.INTER) GO TO 40 IF (IB.EQ.NBNDPT) GO TO 40 C DO NOT LOOK FOR INTERIOR POINT IF THE SUCCESSOR OR C PREDESSOR OF THE BOUNDARY POINT IS ON THE SAME GRID C LINE UNLESS THIS IS A GRID POINT( TYPE = BOTH ) TESTXM = XBOUND(IB) - XBOUND(IB-1) TESTYM = YBOUND(IB) - YBOUND(IB-1) TESTX = XBOUND(IB+1) - XBOUND(IB) TESTY = YBOUND(IB+1) - YBOUND(IB) IF (AMIN1(ABS(TESTX),ABS(TESTXM)).LE.EPSGRD .AND. * TYPE2.EQ.VERT) GO TO 40 IF (AMIN1(ABS(TESTYM),ABS(TESTY)).LE.EPSGRD .AND. * TYPE2.EQ.HORZ) GO TO 40 C DO NOT LOOK FOR INTERIOR POINT IF BOUNDARY IS TANGENT C TO HORIZONTAL OR VERTICAL GRID LINE IF (SIGN(1.,TESTX*TESTXM).LT.0.0 .AND. TYPE2.NE.VERT) GO TO 40 IF (SIGN(1.,TESTY*TESTYM).LT.0.0 .AND. TYPE2.NE.HORZ) GO TO 40 C HAVE FOUND A POINT WHOSE NEIGHBOR SHOULD BE C AN INTERIOR POINT, EXIT LOOP GO TO 70 C UPDATE TYPES FOR NEXT POINT 40 CONTINUE IF (LEVEL.GE.5) WRITE (MOUTPT,99997) IB, TYPE1, TYPE2, TYPE3 99997 FORMAT (5X, 10HTEST POINT, I3, 22H FOR INTERIOR, TYPES =, 3(2X, * A4)) 50 CONTINUE C C DO NOT GIVE UP UNTIL A SECOND PASS IS MADE IF (FIRST) GO TO 20 C *********************** FATAL ERROR ********************* C IF THE LOOP IS ENDED NORMALLY 60 CONTINUE C IF (LEVEL.GE.0) WRITE (MOUTPT,99996) 99996 FORMAT (/5(3H **), 39H FATAL ERROR, NO INITIAL INTERIOR POINT, * /9X, 31HLOCATED NEXT TO BOUNDARY POINTS) FATAL = .TRUE. RETURN C C ---------- NEIGHBOR OF IB POINT SHOULD BE AN INTERIOR POINT --------- C 70 CONTINUE IX = MOD(BGRID(IB),IPACKB) JY = BGRID(IB)/IPACKB C DEBUG OUTPUT FOR LEVELS 3 AND 5 IF (LEVEL.GE.5) WRITE (MOUTPT,99997) IB, TYPE1, TYPE2, TYPE3 IF (LEVEL.GE.3) WRITE (MOUTPT,99995) IB, IX, JY, GTYPE(IX,JY) 99995 FORMAT (9X, 33HCHECK NEIGHBOR OF BOUNDARY POINT , I3, 9H AS INTER, * 17HIOR POINT. IT IS , 2I4, 13H WITH GTYPE , I7) C C TO LOCATE INTERIOR POINT FROM BOUNDARY WE COMPUTE INCREMENTS C IADDX, JADDY FOR (IX,JY) = THE CORNER OF THE GRID SQUARE C CONTAINING THE IB POINT. THESE INCREMENTS ARE DETERMINED C BY THE FOLLOWING( THE SING ORIENT IS USED TO COMBINE THE TWO C ORIENTATIONS) C C CLOCKWISE COUNTERCLOCKWISE C TESTX TESTY TYPE2 IADDX JADDY IADDX JADDY C + HORZ 1 0 C - HORZ 0 1 C + BOTH 1 -1 C - BOTH -1 1 C + BOTH -1 1 C - BOTH 1 -1 C + VERT 1 0 C - VERT 0 1 C C SELECT CASE OF TYPE2 = BOUNDARY POINT TYPE IF (TYPE2.EQ.HORZ) GO TO 90 IF (TYPE2.EQ.VERT) GO TO 100 C HAVE BOUNDARY POINT OF TYPE BOTH, TREAT AS HORZ OR VERT C IF BOUNDARY IS MOVING STEADILY IN X OR Y DIRECTION C CHOOSE THE DIRECTION WITH THE MOST MOVEMENT XMOVE = ABS(TESTX) + ABS(TESTXM) YMOVE = ABS(TESTY) + ABS(TESTYM) XRATIO = AMIN1(ABS(TESTX),ABS(TESTXM))/(XMOVE+EPSGRD) YRATIO = AMIN1(ABS(TESTY),ABS(TESTYM))/(YMOVE+EPSGRD) C DO NOT CHOOSE DIRECTION WHERE MOVEMENT CHANGES SIGN IF (ABS(TESTX+TESTXM).LT.XMOVE-EPSGRD) XRATIO = -1. IF (ABS(TESTY+TESTYM).LT.YMOVE-EPSGRD) YRATIO = -1. RATIO = AMAX1(XRATIO,YRATIO) IF (LEVEL.GE.5) WRITE (MOUTPT,99994) XMOVE, XRATIO, YMOVE, YRATIO 99994 FORMAT (9X, 46HVARIABLES TO CLASSIFY BOTH POINT. XMOVE,XRATIO, * 15H,YMOVE,YRATIO =/15X, 4F12.6) IF (RATIO.GT..01) GO TO 80 C SOMETHING MIGHT BE WRONG, GO BACK TO INITIAL SEARCH GO TO 10 80 IF (RATIO.LE.XRATIO-EPSGRD) GO TO 100 C C IB BOUNDARY POINT IS ON HORIZONTAL LINE 90 IADDX = 0 IF (TESTY*ORIENT.GT.0.0) IADDX = 1 C ADJUST IADDX FOR A POINT OF TYPE BOTH IF (TYPE2.EQ.BOTH .AND. IADDX.EQ.0) IADDX = -1 IX = IX + IADDX C C CHECK THAT TENTATIVE INTERIOR POINT HAS GTYPE THAT POINTS C TO THE IB BOUNDARY POINT IABSG = IABS(GTYPE(IX,JY))/IPACKB IF (LEVEL.GE.5) WRITE (MOUTPT,99993) TYPE2, IADDX, TESTY, ORIENT, * IX, IABSG 99993 FORMAT (9X, A4, I4, 2F12.6, 2I5) C CHECK FOR BOUNDARY TO THE LEFT IF (IADDX.LE.1) IDIFF = MOD(IABSG/8,2) - 1 C CHECK FOR BOUNDARY TO THE RIGHT IF (IADDX.EQ.0) IDIFF = MOD(IABSG/2,2) - 1 IF (IDIFF.EQ.0) GO TO 110 C SOMETHING IS WRONG, GO BACK TO INITIAL SEARCH GO TO 10 C C IB POINT IS ON VERTICAL GRID LINE 100 JADDY = 0 IF (TESTX*ORIENT.LT.0.0) JADDY = 1 IF (TYPE2.EQ.BOTH .AND. JADDY.EQ.0) JADDY = -1 JY = JY + JADDY C C CHECK THAT TENTATIVE INTERIOR POINT HAS GTYPE THAT POINTS C TO THE IB BOUNDARY POINT IABSG = IABS(GTYPE(IX,JY))/IPACKB IF (LEVEL.GE.5) WRITE (MOUTPT,99993) TYPE2, JADDY, TESTX, ORIENT, * JY, IABSG C CHECK FOR BOUNDARY TO THE BELOW IF (JADDY.EQ.1) IDIFF = MOD(IABSG/4,2) - 1 C CHECK FOR BOUNDARY TO THE ABOVE IF (JADDY.LE.0) IDIFF = MOD(IABSG,2) - 1 IF (IDIFF.EQ.0) GO TO 110 C SOMETHING IS WRONG, GO BACK TO INITIAL SEARCH GO TO 10 C C GOT AN INTERIOR POINT 110 CONTINUE IF (LEVEL.GE.3) WRITE (MOUTPT,99992) IX, JY 99992 FORMAT (9X, 10H-----POINT, 2I4, 30H PASSES CHECK TO BE INTERIOR P, * 2HT.) IF (LEVEL.GE.4) WRITE (MOUTPT,99991) IADDX, TESTX, JADDY, TESTY 99991 FORMAT (9X, 30HWITH IADDX,TESTX,JADDY,TESTY =, 2(I3, F12.4)) QX(1) = IX QY(1) = JY GTYPE(IX,JY) = -GTYPE(IX,JY) LENGTH = 1 C ISTART = 1 ISTOP = 1 C PROCESS THE QUEUE BY EXPANDING THE INTERIOR UNTIL NO C POINTS ARE LEFT ON ITS EDGE 120 CONTINUE CALL EXPAND(QX(ISTART), QY(ISTART), NUPTS, INEW, JNEW, GTYPE, * NGDIMX, NGDIMY) C C REMOVE THIS POINT FROM THE QUEUE LENGTH = LENGTH - 1 ISTART = ISTART + 1 IF (ISTART.GT.QLIMIT) ISTART = 1 IF (NUPTS.GT.0) GO TO 130 C C TEST FOR EMPTY QUEUE AND TERMINATION IF (LENGTH.GT.0) GO TO 120 C C *********** SUBPROGRAM EXIT ********** RETURN C 130 DO 140 K=1,NUPTS ISTOP = ISTOP + 1 IF (ISTOP.GT.QLIMIT) ISTOP = 1 LENGTH = LENGTH + 1 C C **** TEST FOR QUEUE OVERFLOW **** IF (LENGTH.GT.QLIMIT) GO TO 150 QX(ISTOP) = INEW(K) QY(ISTOP) = JNEW(K) 140 CONTINUE GO TO 120 C C *********** FATAL ERROR EXIT ****** 150 CONTINUE C IF (LEVEL.GE.0) WRITE (MOUTPT,99990) QLIMIT 99990 FORMAT (/, 5(3H **), 41H FATAL ERROR IN TAGGING ALL INTERIOR POIN, * 2HTS, /, 9X, 31H QUEUE OVERFLOWED IN SUBROUTINE, 12H FILL OF DOM, * 6HAIN , /, 9X, 37H DOMAIN PROCESSOR MUST BE RECOMPILED , * 11HWITH QLIMIT, I5, 10H INCREASED) IF (LEVEL.GE.2) WRITE (MOUTPT,99989) 99989 FORMAT (25X, 23HFINAL QUEUE, THEN ISTOP) IF (LEVEL.GE.2) WRITE (MOUTPT,99988) (K,QX(K),QY(K),K=1,QLIMIT), * ISTOP 99988 FORMAT (6(4X, 3I3)) FATAL = .TRUE. RETURN END SUBROUTINE GVALUS(XGRID, YGRID, NGDIMX, NGDIMY, GTYPE, NBDIM, GVA 10 * XBOUND, YBOUND, BPTYPE, BGRID) C C *** THIS PROGRAM LOCATES THE NEIGHBORING GRID POINTS WITH RESPECT C TO THE BOUNDARY. THE INFO IS ENCODED INTO GTYPE. C GTYPE IS INITIALLY NEGATIVE BUT IT IS SET POSITIVE AT C INTERIOR POINTS BY THE SUBROUTINE FILL CALLED LATER. C C INPUT. XBOUND,YBOUND = BOUNDARY POINT COORDINATES C BGRID = - - GRID LOCATIONS C BPTYPE = - - TYPES C XGRID,YGRID = GRID DEFINITION C GTYPE = GRID MARKERS INITIALIZED TO ZERO C OUTPUT. GTYPE MODIFIED C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C REAL XGRID(NGDIMX), YGRID(NGDIMY) DIMENSION GTYPE(NGDIMX,NGDIMY) DIMENSION XBOUND(NBDIM), YBOUND(NBDIM), BPTYPE(NBDIM), * BGRID(NBDIM) C C LOOP OVER BOUNDARY POINTS DO 120 NB=N1BDPT,NBNDPT IX = MOD(BGRID(NB),IPACKB) JY = BGRID(NB)/IPACKB C NBPTDL = NBNDPT - N1BDPT + 1 NBACK = MOD(NB-N1BDPT+1+NBPTDL-2,NBPTDL) + N1BDPT C ELIMINATE ANY INTERIOR POINT BEHIND THE NB POINT IF (BPTYPE(NBACK).EQ.INTER) NBACK = MOD(NB-N1BDPT+1+NBPTDL-3, * NBPTDL) + N1BDPT C VARIABLES FOR TESTING THE DIRECTION OF THE BOUNDARY TESTX = XBOUND(NB) - XBOUND(NBACK) TESTY = YBOUND(NB) - YBOUND(NBACK) TYPE = BPTYPE(NB) C C THERE ARE SEVERAL CASES = VERT, HORZ AND BOTH,INTER IF (TYPE.EQ.VERT) GO TO 10 IF (TYPE.EQ.HORZ) GO TO 30 IF (TYPE.EQ.BOTH) GO TO 50 IF (TYPE.EQ.INTER) GO TO 100 C C HAVE ERROR IN ILLEGAL TYPE C IF (LEVEL.GE.0) WRITE (MOUTPT,99999) NB, IX, JY, TYPE 99999 FORMAT (/5(3H **), 40H FATAL ERROR IN DETERMINING NEIGHBORS OF, * /9X, 15HBOUNDARY POINT , I3, 16H IN GRID SQUARE , 2I3, /9X, * 21H WITH UNKNOWN TYPE = , A4) FATAL = .TRUE. GTYPE(IX,JY) = NB GO TO 110 C C ****** BOUNDARY POINT IS ON VERTICAL GRID LINE ************** 10 CONTINUE C CHECK FOR A DOUBLE CROSSING ON THIS GRID IF (ABS(TESTX).LE.EPSGRD .AND. YBOUND(NBACK).LT.YGRID(JY+1) * +EPSGRD .AND. YBOUND(NBACK).GT.YGRID(JY)-EPSGRD) GO TO 20 C C SET GTYPE FOR POINTS ABOVE AND BELOW GTYPE(IX,JY) = ISETGT(GTYPE(IX,JY),1,NB) GTYPE(IX,JY+1) = ISETGT(GTYPE(IX,JY+1),3,NB) GO TO 110 C C HAVE A DOUBLE CROSSING C GTYPE HAS ALREADY BEEN SET CORRECTLY BY THE NBACK POINT C EXCEPT AT THE START OF THE BOUNDARY 20 IF (NB.GT.2) GO TO 110 C C IDIREC = +1 FOR BOUNDARY GOING UP, = 0 GOING DOWN IDIREC = +1 IF (TESTY.LT.0.0) IDIREC = 0 JY1ID = JY + 1 - IDIREC GTYPE(IX,JY1ID) = -NBACK - IPACKB*(4-3*IDIREC) GO TO 110 C C ******* BOUNDARY POINT IS ON HORIZONTAL GRID LINE *********** 30 CONTINUE C CHECK FOR A DOUBLE CROSSING ON THIS GRID IF (ABS(TESTY).LE.EPSGRD .AND. XBOUND(NBACK).LT.XGRID(IX+1) * +EPSGRD .AND. XBOUND(NBACK).GT.XGRID(IX)-EPSGRD) GO TO 40 C C SET GTYPE FOR POINTS LEFT AND RIGHT GTYPE(IX,JY) = ISETGT(GTYPE(IX,JY),2,NB) GTYPE(IX+1,JY) = ISETGT(GTYPE(IX+1,JY),4,NB) GO TO 110 C C HAVE A DOUBLE CROSSING C GTYPE HAS ALREADY BEEN SET CORRECTLY BY THE NBACK POINT C EXCEPT AT THE START OF THE BOUNDARY 40 IF (NB.GT.2) GO TO 110 C C IDIREC = +1 FOR BOUNDARY GOING RIGHT, = 0 GOING LEFT IDIREC = +1 IF (TESTX.LT.0.0) IDIREC = 0 IX1ID = IX + 1 - IDIREC GTYPE(IX1ID,JY) = -NBACK - IPACKB*(2+6*IDIREC) GO TO 110 C C ******* BOUNDARY POINT IS AT GRID POINT IX,JY ***************** 50 CONTINUE C C TYPE THIS GRID POINT AS A BOUNDARY POINT GTYPE(IX,JY) = NB C SET ALL 4 NEIGHBORS THAT ARE NOT ALSO BOUNDARY POINTS OF C THE GRID. SKIP THOSE BOUNDARY PTS THAT ARE GRID PTS OR C HAVE BOUNDARY PTS ALREADY INDICATED ON THIS SEGMENT IF (IX.EQ.NGRIDX) GO TO 60 C POINT TO THE RIGHT GTYPE(IX+1,JY) = ISETGT(GTYPE(IX+1,JY),4,NB) 60 IF (JY.EQ.NGRIDY) GO TO 70 C POINT ABOVE GTYPE(IX,JY+1) = ISETGT(GTYPE(IX,JY+1),3,NB) 70 IF (JY.EQ.1) GO TO 80 C POINT BELOW GTYPE(IX,JY-1) = ISETGT(GTYPE(IX,JY-1),1,NB) 80 IF (IX.EQ.1) GO TO 90 C POINT TO THE LEFT GTYPE(IX-1,JY) = ISETGT(GTYPE(IX-1,JY),2,NB) 90 CONTINUE GO TO 110 C C BOUNDARY POINT IS INTERIOR, NO GTYPES TO SET 100 CONTINUE C C ALL CASES COME HERE, POSSIBLE PRINT FOR DEBUG 110 CONTINUE C C LEVEL 4 DEBUG OUTPUT C NOTE: GTYPE VALUES MAY BE PRINTED FOR ILLEGAL INDICES C 0, NGRIDX + 1, NGRIDY + 1 IF (LEVEL.GE.4) WRITE (MOUTPT,99998) NB, TYPE, GTYPE(IX,JY), * GTYPE(IX+1,JY), GTYPE(IX-1,JY), GTYPE(IX,JY+1), GTYPE(IX,JY-1) 99998 FORMAT (5X, 5HPOINT, I4, 9H OF TYPE , A4, 15H AND GTYPE = , * 5I7) IF (LEVEL.GE.5) WRITE (MOUTPT,99997) TESTX, TESTY, NBPTDL, NBACK 99997 FORMAT (14X, 27HTESTX,TESTY,NBDPDL,NBACK = , 2F12.7, I4, I7) 120 CONTINUE C END OF LOOP OVER BOUNDARY C IF (LEVEL.LT.2) RETURN WRITE (MOUTPT,99996) 99996 FORMAT (9X, 34H----HAVE SET GTYPE FOR ALL POINTS ) RETURN END SUBROUTINE HOLE(XGRID, YGRID, NGDIMX, NGDIMY, BRANGE, NPDIM, HOL 10 * BCOORD, SCLOCK, SARC, SLEVEL, GTYPE, XBOUND, YBOUND, PIECE, * BPTYPE, BNEIGH, BGRID, BPARAM, NBDIM, NBPTS, FAIL, GTYPE2) C C ***** THIS ROUTINE IS AN ALTERNATE TO REGION AS A DRIVER FOR THE C DOMAIN PROCESSOR THAT IS USED TO REMOVE HOLES FROM THE DOMAIN. C THE HOLE IS PROCESSED LIKE AN ORIGINAL DOMAIN AND THEN THE C RESULTING INFORMATION USED BY REMOVH TO REMOVE THE HOLE. C THE INPUT/OUTPUT IS THE SAME AS REGION - SEE COMMENTS THERE - C EXCEPT THAT A SECOND ARRAY( = GTYPE2 ) IS PASSED AS C WORKSPACE FOR THE NEW GRID POINT TYPES C ***** REGION MUST BE CALLED BEFORE HOLE ***** C REAL XGRID(NGDIMX), YGRID(NGDIMY), BRANGE(2,NPDIM) LOGICAL SCLOCK, SARC, FAIL INTEGER SLEVEL, GTYPE2(NGDIMX,NGDIMY) DIMENSION GTYPE(NGDIMX,NGDIMY), XBOUND(NBDIM), YBOUND(NBDIM), * BPTYPE(NBDIM), BNEIGH(NBDIM), BGRID(NBDIM), BPARAM(NBDIM), * PIECE(NBDIM) EXTERNAL BCOORD INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C IF (LEVEL.GT.0) WRITE (MOUTPT,99999) 99999 FORMAT (///1H , 19(1H-)/1H , 18H DOMAIN PROCESSOR/1H , 19(1H-)// * 5X, 29H H O L E P R O C E S S O R ) C C RESET COMMON BLOCK VARIABLES THAT MAY HAVE BEEN CHANGED C FROM THE PRECEEDING CALL TO REGION OR HOLE C NBOUND = NPDIM LEVEL = SLEVEL CLOCKW = SCLOCK ARC = SARC C C INITIALIZE JUMP TO HOLE, BOUNDARY POINT COUNT BPTYPE(NBNDPT+1) = JUMP N1BDPT = NBNDPT + 2 C MARK THAT WE ARE IN HOLE INHOLE = .TRUE. C CALL DOMAIN(XGRID, YGRID, NGDIMX, NGDIMY, BRANGE, NPDIM, BCOORD, * GTYPE2, XBOUND, YBOUND, PIECE, BPTYPE, BNEIGH, BGRID, BPARAM, * NBDIM) C FAIL = FATAL IF (FATAL) RETURN C C USE GTYPE2 TO REMOVE HOLE FROM DOMAIN CALL REMOVH(GTYPE, GTYPE2, NGDIMX, NGDIMY) FAIL = FATAL N1BND = NBOUND + 1 NBPTS = NBNDPT C RETURN END LOGICAL FUNCTION INSIDE(X, Y, IX, JY, INEXT, JNEXT, XGRID, INS 10 * NGDIMX, YGRID, NGDIMY) C C **** THIS FUNCTION TESTS WHETHER THE POINT (X,Y) IS INSIDE THE GRID C WITH SIDES IX AND INEXT, JY AND JNEXT C THIS HAS BEEN MADE INTO A FUNCTION SUBPROGRAM BECAUSE TOO C MANY COMPILERS FAILED ON IT AS AN ARITHMETIC STATEMENT C FUNCTION. C REAL X, Y, XGRID(NGDIMX), YGRID(NGDIMY) INTEGER IX, JY, INEXT, JNEXT C C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C INSIDE = X.GE.(XGRID(IX)-EPSGRD) .AND. X.LE.(XGRID(INEXT)+EPSGRD) * .AND. Y.GE.(YGRID(JY)-EPSGRD) .AND. Y.LE.(YGRID(JNEXT)+EPSGRD) RETURN END INTEGER FUNCTION ISETGT(GTYPE, IPOINT, NB) ISE 10 C C **** THIS FUNCTION SETS THE GTYPE VALUE OF A GRID POINT. C 1. SETS IT = INSIDE FOR INTERIOR POINTS C 2. UPDATES IT IF GTYPE HAS BEEN SET BY SOME OTHER BOUNDARY PT C 3. DOES NOTHING IF GRID POINT IS A BOUNDARY POINT C 4. SETS ITS FIRST VALUE IF IT HAS NOT BEEN SET BEFORE C LOGICAL PTNAY INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C C PTNAY IS TRUE FOR BOUNDARY POINT LOCATED C ABOVE J=1 C RIGHT J=2 C BELOW J=3 C LEFT J=4 PTNAY(K,J) = MOD(K/(IPACKB*2**(MOD(J-1,4))),2).EQ.1 C C ISETGT = GTYPE C DO NOT CHANGE GTYPE AT BOUNDARY POINTS OR AT DOUBLE CROSSINGS IF ((ISETGT.LT.NSIDE .AND. ISETGT.GT.0) .OR. PTNAY(IABS(ISETGT), * IPOINT)) RETURN ISETGT = IABS(GTYPE) IF (ISETGT.EQ.0) ISETGT = NSIDE ISETGT = -MIN0(NB,MOD(ISETGT,IPACKB)) - IPACKB*(ISETGT/IPACKB+2** * (IPOINT-1)) RETURN END SUBROUTINE LOCATE(XB, YB, IX, JY, TYPEC, XGRID, YGRID, NGDIMX, LOC 10 * NGDIMY) C C *** THIS PROGRAM LOCATES AND TYPES A GIVEN BOUNDARY POINT XB,YB C IT FINDS THE GRID SQUARE BOUNDED BY X-GRID LINES IX,IX+1 AND C Y-GRID LINES JY,JY+1 WHICH CONTAIN TH XB,YB POOINT. C INPUT = XB,YB = POINT COORDINATES C XGRID,YGRID = GRID DEFINITION C OUTPUT= IX,JY = GRID CONTAINING XB,YB C TYPEC = TYPE, BUT NOT RETURNED THRU COMMON C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C INTEGER TYPEC, TYPE1, TYPE2 REAL XGRID(NGDIMX), YGRID(NGDIMY) C C FIND X-GRID LINE TO LEFT OF XB IX = 0 DO 10 IXX=1,NGRIDX IX = IX + 1 IF (XB.LE.XGRID(IX)-EPSGRD) GO TO 20 10 CONTINUE IF (XB.LE.XGRID(NGRIDX)+EPSGRD) GO TO 20 C IF (LEVEL.GE.0) WRITE (MOUTPT,99999) XB, XGRID(NGRIDX) 99999 FORMAT (/5(3H **), 41H FATAL ERROR, BOUNDARY POINT X-COORDINATE/ * 9X, F10.6, 20H BIGGER THAN X-GRID , F10.6) FATAL = .TRUE. RETURN 20 IX = IX - 1 IF (XB.GT.XGRID(1)-EPSGRD) GO TO 30 C WRITE (MOUTPT,99998) XB, XGRID(1) 99998 FORMAT (/5(3H **), 41H FATAL ERROR, BOUNDARY POINT X-COORDINATE/ * 9X, F10.6, 18H LESS THAN X-GRID , F10.6) FATAL = .TRUE. RETURN C C C FIND Y-GRID LINE BELOW YB 30 CONTINUE JY = 0 DO 40 JYY=1,NGRIDY JY = JY + 1 IF (YB.LE.YGRID(JY)-EPSGRD) GO TO 50 40 CONTINUE IF (YB.LE.YGRID(NGRIDY)+EPSGRD) GO TO 50 C WRITE (MOUTPT,99997) YB, YGRID(NGRIDY) 99997 FORMAT (/5(3H **), 41H FATAL ERROR, BOUNDARY POINT Y-COORDINATE/ * 9X, F10.6, 20H BIGGER THAN Y-GRID , F10.6) FATAL = .TRUE. RETURN 50 JY = JY - 1 IF (YB.GT.YGRID(1)-EPSGRD) GO TO 60 C WRITE (MOUTPT,99996) YB, YGRID(1) 99996 FORMAT (/5(3H **), 41H FATAL ERROR, BOUNDARY POINT Y-COORDINATE/ * 9X, F10.6, 18H LESS THAN Y-GRID , F10.6) FATAL = .TRUE. RETURN C C DETERMINE IF FIRST POINT IS ON GRID LINES 60 TYPE1 = INTER C CHECK X-GRID LINES IF (ABS(XB-XGRID(IX)).LE.EPSGRD) TYPE1 = VERT C CHECK NEXT X-GRID LINE IF (ABS(XB-XGRID(IX+1)).GT.EPSGRD) GO TO 70 TYPE1 = VERT IX = IX + 1 C CHECK Y-GRID LINES 70 TYPE2 = INTER IF (ABS(YB-YGRID(JY)).LE.EPSGRD) TYPE2 = HORZ C CHECK NEXT Y-GRID LINE IF (ABS(YB-YGRID(JY+1)).GT.EPSGRD) GO TO 80 TYPE2 = HORZ JY = JY + 1 C CHECK FOR BOTH 80 CONTINUE TYPEC = INTER IF (TYPE1.EQ.INTER .AND. TYPE2.EQ.HORZ) TYPEC = HORZ IF (TYPE1.EQ.VERT .AND. TYPE2.EQ.INTER) TYPEC = VERT IF (TYPE1.EQ.VERT .AND. TYPE2.EQ.HORZ) TYPEC = BOTH C DEBUG IF (LEVEL.GE.4) WRITE (MOUTPT,99995) XB, YB, IX, JY, TYPEC, * TYPE1, TYPE2 99995 FORMAT (5X, 14HLOCATED POINT , 2F10.6, 9H IN GRID , 2I3, 6H OF TY, * 3HPE , A4, 14H AUX. VARS. = , A4, 1X, A4) RETURN END SUBROUTINE NEIGH(XGRID, YGRID, NGDIMX, NGDIMY, GTYPE, NBDIM, NEI 10 * XBOUND, YBOUND, BNEIGH, BPTYPE, BGRID) C C ***** THIS ROUTINE COMPUTES POINTERS FROM THE BOUNDARY POINTS TO C THE INTERIOR GRID POINTS. THE VALUE IS ENCODED INTO BNEIGH C THE SAME WAY POINTERS ARE ENCODED INTO GTYPE( SEE COMMENTS IN C DOMAIN ) C C INPUT. XGRID,YGRID = GRID DEFINITION C XBOUND,YBOUND = BOUNDARY POINT COORDINATES C BPTYPE = - - TYPES C BGRID = - - GRID LOCATIONS C GTYPE = GRID POINT MARKERS C OUTPUT. BNEIGH = POINTERS FROM BOUNDARY TO GRID C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C REAL XGRID(NGDIMX), YGRID(NGDIMY) DIMENSION GTYPE(NGDIMX,NGDIMY) DIMENSION XBOUND(NBDIM), YBOUND(NBDIM), BPTYPE(NBDIM), * BGRID(NBDIM), BNEIGH(NBDIM) C IF CALLED FROM HOLE WE ANTICIPATE THE CHANGE IN SIGN C OF GTYPE THAT IS COMING LATER BY SETTING IS = -1 IS = 1 IF (INHOLE) IS = -1 C INITIALIZE POINTER TO PREVIOUS JUMP POINT LASTJ = N1BDPT - 1 C C LOOP OVER BOUNDARY SETTING BNEIGH VALUES DO 90 IB=N1BDPT,NBNDPT C FIND GRID CONTAINING BOUNDARY POINT IB IX = MOD(BGRID(IB),IPACKB) JY = BGRID(IB)/IPACKB C C SELECT CASE DEPENDING ON TYPE OF POINT TYPE = BPTYPE(IB) BNEIGH(IB) = 0 IF (TYPE.EQ.HORZ) GO TO 10 IF (TYPE.EQ.VERT) GO TO 20 IF (TYPE.EQ.BOTH) GO TO 30 IF (TYPE.EQ.JUMP) GO TO 70 C EVERYTHING ELSE HAS BNEIGH = 0 GO TO 80 C C HORIZONTAL GRID POINT 10 CONTINUE IF (IS*GTYPE(IX,JY).GE.NSIDE) BNEIGH(IB) = BNEIGH(IB) + 8 IF (IX.EQ.NGRIDX) GO TO 80 IF (IS*GTYPE(IX+1,JY).GE.NSIDE) BNEIGH(IB) = BNEIGH(IB) + 2 GO TO 80 C C VERTICAL GRID POINT 20 CONTINUE IF (IS*GTYPE(IX,JY).GE.NSIDE) BNEIGH(IB) = BNEIGH(IB) + 4 IF (JY.EQ.NGRIDY) GO TO 80 IF (IS*GTYPE(IX,JY+1).GE.NSIDE) BNEIGH(IB) = BNEIGH(IB) + 1 GO TO 80 C C BOTH GRID POINT 30 CONTINUE IF (JY.EQ.NGRIDY) GO TO 40 IF (IS*GTYPE(IX,JY+1).GE.NSIDE) BNEIGH(IB) = BNEIGH(IB) + 1 40 IF (IX.EQ.NGRIDX) GO TO 50 IF (IS*GTYPE(IX+1,JY).GE.NSIDE) BNEIGH(IB) = BNEIGH(IB) + 2 50 IF (JY.EQ.1) GO TO 60 IF (IS*GTYPE(IX,JY-1).GE.NSIDE) BNEIGH(IB) = BNEIGH(IB) + 4 60 IF (IX.EQ.1) GO TO 80 IF (IS*GTYPE(IX-1,JY).GE.NSIDE) BNEIGH(IB) = BNEIGH(IB) + 8 GO TO 80 C JUMP POINT, COPY I1BNGH AND RESET LASTJ 70 BNEIGH(IB) = BNEIGH(LASTJ+1) LASTJ = IB C C TERMINAL MESSAGE 80 IF (LEVEL.GE.4) WRITE (MOUTPT,99999) IB, IX, JY, TYPE, * BNEIGH(IB) 90 CONTINUE 99999 FORMAT (5X, 13HNEIGHBORS OF , I3, 11H IN SQUARE , 2I3, 8H OF TYPE, * 1H , A4, 16H FOUND, BNEIGH =, I7) C C PROGRAM EXIT C SET BNEIGH FOR LAST+1 POINT BNEIGH(NBNDPT+1) = BNEIGH(LASTJ+1) IF (LEVEL.GE.3) WRITE (MOUTPT,99998) 99998 FORMAT (9X, 32HHAVE SET BNEIGH FOR THE BOUNDARY) RETURN END SUBROUTINE REGULA(TARGET, XCASE, GESSL, GESSR, PSTOP, EPSILN, REG 10 * BCOORD) C C ***** THIS IS A SIMPLE MODIFIED REGULA FALSI METHOD TO IMPROVE THE C GUESSES FOR SECANT. THIS SPECIAL REGULA FALSI USES THE USUAL C FORMULA IS THE NEXT POINT IS BETWEEN THE PREVIOUS PAIR AND C USES A PROPORTIONAL RULE OTHERWISE WHICH KEEPS GESSL FIXED C IT IS NOT RUN TO COMPLETE CONVERGENCE IN NORMAL USE, THE SECANT C METHOD IS USED FOR THAT. C EPSILN = CONVERGENCE TOLERANCE TIMES 10. C CONVERG = EPSILN*.1 C ITERATIONS LIMITED TO 6 WHEN TARGET NOT BRACKETED C LIMITED TO 20 OTHERWISE C GESSL = FIRST GUESS OF INCREMENT TO TARGET C START = CORRESPONDING VARIABLE VALUE - NEVER CHANGED C VLEFT = CORRESPONDING VARIABLE VALUE C ELEFT = CORRESPONDING ERROR C GESSR = SECOND GUESS OF INCREMENT TO TARGET C VRIGHT= CORRESPONDING VARIABLE VALUE C ERIGHT= CORRESPONDING ERROR C GESNEW = NEW VALUE OF GESSL FROM REGULA FALSI FORMULA C VNEW = CORRESPONDING VARIABLE VALUE C ENEW = CORRESPONDING ERROR C DELJIG = SMALL VALUE TO JIGGLE RESULTS TO AVOID DUPLICATES C PSTOP = LIMIT ON THE PARAMETER VALUE C DELPMX = LIMIT ON THE CHANGE FROM PARAM = PSTOP-PARAM C GTOP = TOP LIMIT ON PARAMETER CHANGE C GBOT = BOTTOM LIMIT ON PARAMETER CHANGE C XCASE = .TRUE. WHEN TARGET IS X-VALUE C .FALSE. WHEN TARGET IS Y-VALUE C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C EXTERNAL BCOORD LOGICAL XCASE C FIND COORDINATES, ERRORS FOR INITIAL GUESSES CALL BCOORD(PARAM+GESSL, XLEFT, YLEFT, IPIECE) START = YLEFT IF (XCASE) START = XLEFT VLEFT = START ELEFT = TARGET - VLEFT CALL BCOORD(PARAM+GESSR, XRIGHT, YRIGHT, IPIECE) VRIGHT = YRIGHT IF (XCASE) VRIGHT = XRIGHT ERIGHT = TARGET - VRIGHT ENEW = 0.0E0 C SET LIMIT FOR PARAMETER DURING THE SEARCH DELPMX = PSTOP - PARAM C SIGN CHANGE BETWEEN GUESSES, RESET DELPMX IF (ELEFT*ERIGHT.LT.0.) DELPMX = GESSR CONVRG = .1E0*EPSILN GBOT = GESSL GTOP = DELPMX C C START MODIFIED REGULA FALSI ITERATION LOOP KREG = 0 10 KREG = KREG + 1 C PROTECTION FROM DIVISION BY ZERO DENOM = VLEFT - VRIGHT IF (ABS(DENOM).LE.EPSGRD) GO TO 40 C REGULA FALSI FORMULA GESNEW = GESSR - ERIGHT*(GESSR-GESSL)/DENOM C CHECK BRACKETEING OF NEW POINT IF (ABS(GESSR-GESSL).LT.ABS(GESSR-GESNEW)+ABS(GESSL-GESNEW)) GO * TO 50 C HAVE GESNEW BRACKETED, PROCEED WITH MODIFIED REGULA FALSI CALL BCOORD(PARAM+GESNEW, XNEW, YNEW, IPIECE) VNEW = YNEW IF (XCASE) VNEW = XNEW EOLD = ENEW ENEW = TARGET - VNEW IF (ENEW*ERIGHT.GT.0.0) GO TO 20 C REPLACE VLEFT WITH NEW POINT VLEFT = VNEW ELEFT = ENEW GESSL = GESNEW IF (ENEW*EOLD.GT.0.) ERIGHT = ERIGHT/2.E0 GO TO 30 C REPLACE VRIGHT POINT WITH NEW POINT 20 VRIGHT = VNEW ERIGHT = ENEW GESSR = GESNEW IF (ENEW*EOLD.GT.0.) ELEFT = ELEFT/2.E0 30 CONTINUE IF (KREG.GE.20 .OR. ABS(ENEW).LE.CONVRG) GO TO 60 GO TO 10 C C ABOUT TO DIVIDE BY ZERO, CONJURE UP NEW DENOM 40 DENOM = TARGET - VLEFT + FLOAT(KREG)*EPSGRD*5. C DO NOT HAVE SOLUTION BRACKETED, DO NOT MOVE C STARTING POINT( GESSL ), ADJUST GESSR 50 GESSR = (GESSL*ELEFT-ERIGHT*GESSR)/DENOM C JIGGLE SOME TO AVOID GETTING DUPLICATE POINTS DELJIG = DELPMX*.01E0/FLOAT(KREG) GESSR = AMIN1(GTOP-DELJIG,AMAX1(GBOT,DELJIG)) CALL BCOORD(PARAM+GESSR, XRIGHT, YRIGHT, IPIECE) VRIGHT = YRIGHT IF (XCASE) VRIGHT = XRIGHT ERIGHT = TARGET - VRIGHT IF (KREG.LT.6) GO TO 10 C C ERROR EXIT - DID NOT GET A SATISFACTORY IMPROVEMENT C SET GESSR = PROPORTIONAL DISTANCE TO PSTOP CALL BCOORD(PSTOP, XRIGHT, YRIGHT, IPIECE) VRIGHT = YRIGHT IF (XCASE) VRIGHT = XRIGHT GESSR = GESSL - (TARGET-VRIGHT)*(VRIGHT-VLEFT)/(DELPMX-GESSL) DELJIG = (PSTOP-START)*.001 GESSR = AMIN1(GTOP-DELJIG,AMAX1(GESSR,DELJIG)) IF (LEVEL.GE.4) WRITE (MOUTPT,99999) PARAM, TARGET, GESSL, GESSR, * VRIGHT 99999 FORMAT (3(3H ++), 32HREGULA DID NOT CONVERGE NORMALLY, 8H, STARTE, * 14HD WITH PARAM =, F10.6/9X, 13HFOR TARGET = , F10.6, 8H ENDED W, * 29HITH GESSL, GESSR AND VALUE = , 3F11.6) RETURN C C NORMAL EXIT 60 CONTINUE IF (LEVEL.GE.4) WRITE (MOUTPT,99998) START, KREG, GESSL, GESSR, * VRIGHT, VLEFT 99998 FORMAT (9X, 45HREGULA-FALSI START,KREG,GESSL,GESSR,VALUES = , * F10.6, I3, 4F10.6) RETURN END SUBROUTINE REMOVH(GTYPE1, GTYPE2, NDX, NDY) REM 10 C INTEGER GTYPE1(NDX,NDY), GTYPE2(NDX,NDY) LOGICAL FATALE C C C *** REMOVH REMOVES A HOLE FROM A DOMAIN. GTYPE1 REPRESENTS THE C DOMAIN AND GTYPE2 THE HOLE. ON EXIT, GTYPE1 HAS THE HOLE C REMOVED FROM IT. C C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C FATALE = .FALSE. DO 20 I=1,NDX DO 10 J=1,NDY IF (GTYPE2(I,J).EQ.0) GO TO 10 IF (GTYPE1(I,J).NE.NSIDE) FATALE = .TRUE. GTYPE1(I,J) = -GTYPE2(I,J) IF (IABS(GTYPE1(I,J)).LT.NSIDE) GTYPE1(I,J) = * IABS(GTYPE1(I,J)) IF (GTYPE1(I,J).EQ.-NSIDE) GTYPE1(I,J) = 0 10 CONTINUE 20 CONTINUE C IF (.NOT.FATALE) GO TO 30 C C IF (LEVEL.GE.0) WRITE (MOUTPT,99999) 99999 FORMAT (/5(3H **), 32H FATAL ERROR IN ROUTINE REMOVH, /9X, * 36HHOLE TOO NEAR BOUNDARY OR OTHER HOLE/9X, 17HMUST ADJUST GRID , * 16HOR MAKE IT FINER) FATAL = .TRUE. C 30 CONTINUE IF (LEVEL.GT.1) CALL TABLGT(MOUTPT, 1, NDX, NDY, GTYPE1) C RETURN END SUBROUTINE SECANT(PFIRST, PSECND, TARGET, XCASE, PANS, X2, Y2, SEC 10 * OUTCOM, PSTOP, VGRID, BCOORD) C C ******* THIS IS A GENERAL PURPOSE SECANT ROUTINE FOR BOUNDARY C INTERSECTIONS. XCASE TELLS WHETHER AN X OR Y INTERSECTION C IS SOUGHT. CONVERGENCE TEST IS .2*EPSGRD BUT POINT ACCEPTED C AT END EVEN WITHIN ONLY EPSGRD OF TARGET. C OUTCOM CAN BE = SUCCES OR FAIL C C INPUT. PFIRST,PSECOND = GUESSES FOR PARAMETER C PSTOP = LIMIT ON PARAMETER C TARGET,XCASE = DESIRED COORDINATE VALUE( XCASE = .TRUE. C FOR X-VALUE, .FALSE. FOR Y-VALUE ) C VGRID = GRID SIZE FOR THE UNKNOWN VARIABLE C BCOORD = FUNCTION THAT DEFINES BOUNDARY C OUTPUT. PANS,X2,Y2 = ANSWERS( PARAMETER AND X,Y VALUES) C OUTCOM = OUTCOME FLAG ( = SUCCES OR FAIL ) C C INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, NGRIDY, * LEVEL, QLIMIT, QX(250), QY(250), GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON /DMCINT/ N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, NGRIDX, * NGRIDY, LEVEL, QLIMIT, QX, QY COMMON /DMCREL/ PARAM, EPSGRD, EPSTAN COMMON /BNDRYI/ IPIECE, NBOUND, NBNDPT COMMON /BNDRYL/ CLOCKW, ARC, FATAL, INHOLE COMMON /DMCHAR/ TYPE, HORZ, VERT, BOTH, INTER, JUMP C EXTERNAL BCOORD LOGICAL XCASE, LAST INTEGER OUTCOM, SUCCES, FAIL DATA SUCCES, FAIL /4HSUCC,4HFAIL/ C OUTCOM = SUCCES C C SIMPLE PROTECTION AGAINST BAD INITIAL GUESSES DELV = PSECND - PFIRST IF (ABS(DELV).GE.EPSGRD) GO TO 10 IF (DELV.EQ.0.0E0) DELV = .0001E0*(PSTOP-PARAM) PSECND = PSECND + 2.E0*DELV PFIRST = PFIRST - 2.E0*DELV 10 CONTINUE P1 = AMAX1(PARAM,AMIN1(PFIRST,PSTOP)) P2 = AMAX1(PARAM,AMIN1(PSECND,PSTOP)) CALL BCOORD(P1, X1, Y1, IPIECE) CALL BCOORD(P2, X2, Y2, IPIECE) C C INTRODUCE VARIABLES V1,V2 AS UNKNOWNS FOR THE SECANT IF (XCASE) GO TO 20 V1 = Y1 V2 = Y2 GO TO 30 20 V1 = X1 V2 = X2 30 CONTINUE C C SET BOUNDS ON RANGE OF UNKNOWN DELV = TARGET + SIGN(1.,TARGET-V1)*VGRID VLOW = AMIN1(V1,DELV) - EPSGRD VHIGH = AMAX1(V1,DELV) + EPSGRD IF (LEVEL.GE.5) WRITE (MOUTPT,99999) VLOW, VHIGH 99999 FORMAT (9X, 26HBOUNDS ON VARIABLE RANGE =, F12.5, 4H TO , F12.5) IF (LEVEL.GE.5) WRITE (MOUTPT,99998) P1, V1, P2, V2, TARGET, XCASE 99998 FORMAT (9X, 16HSTART SECANT AT , 2(F14.7, F12.7), 12H FOR TARGET , * F11.6, 14H WITH XCASE = , L3) C C SET CONVERGENCE TOLERANCE LESS THAN EPSGRD C THE FACTOR 0.2 MAY BE CHANGED TO TUNE THE CONVERGENCE CONVRG = 0.2E0*EPSGRD ERPASS = EPSGRD EPS01 = EPSGRD*.01E0 C LIMIT THE SECANT METHOD TO 12 ITERATIONS C THERE IS A SPECIAL LAST ITERATION LAST = .FALSE. DO 70 KNTSEC=1,12 C PROTECT AGAINST DIVISION BY ZERO IN SECANT FORMULA IF (ABS(V2-V1).GT.EPS01) PANS = P2 - (V2-TARGET)*(P1-P2)/(V1-V2) IF (ABS(V2-V1).LE.EPS01) PANS = P2 - (V2-TARGET)*(P1-P2)/ * (EPSGRD*SIGN(1.E0,V1-V2)) C RESTRICT PARAMETER TO VALID RANGE C INTRODUCE TEMPORARY PARAMETER VALUE PTEMP PTEMP = AMIN1(PSTOP,AMAX1(PANS,PARAM)) CALL BCOORD(PTEMP, XNU, YNU, IPIECE) VNU = YNU IF (XCASE) VNU = XNU C C CHECK THAT THE SECANT METHOD IS STAYING WITHIN RANGE IF (VNU.LE.VHIGH .AND. VNU.GE.VLOW) GO TO 40 C CHANGE PANS TO KEEP IT NEAR THE LIMIT POINT DELV = VNU - V2 VF = VHIGH - V2 IF (VNU.LT.VLOW) VF = VLOW - V2 RATIO = .4 + .4/FLOAT(KNTSEC) IF (ABS(DELV).GT.EPS01) RATIO = VF/(DELV) PANS = RATIO*(PANS-P2) + P2 IF (LEVEL.GE.5) WRITE (MOUTPT,99997) KNTSEC, P1, P2, PANS, V1, * V2, VNU 99997 FORMAT (15X, 39HRANGE LIMIT INFO P1,P2,PANS,V1,V2,VNU =, I3, * 6F12.7) C RESTRICT PANS TO VALID RANGE C JIGGLED BY .02/KNTSEC TO PREVENT GETTING C DUPLICATE POINTS VF = 1. - .02E0/FLOAT(KNTSEC) PANS = AMIN1(PSTOP*VF,AMAX1(PANS,PARAM)/VF) CALL BCOORD(PANS, XNU, YNU, IPIECE) VNU = YNU IF (XCASE) VNU = XNU GO TO 50 C OK TO USE PTEMP AS PANS 40 PANS = PTEMP C C GET CONVERGENCE CLOSER THAN EPSGRD 50 ERROR = ABS(VNU-TARGET) IF (ERROR.GE.ERPASS) GO TO 60 C SAVE BEST POINT CLOSER THAN EPSGRD ERPASS = ERROR XPASS = XNU YPASS = YNU IF (ERROR.GT.CONVRG) GO TO 60 C HAVE PASSED THE CONVERGENCE TEST C CHECK FOR THE LAST ITERATION IF (LAST) GO TO 90 LAST = .TRUE. C C AFTER NORMAL CONVERGENCE DO 1 MORE ITERATION BEFORE QUITTING C WILL PICK THE BEST OF THE TWO RESULTS CONVRG = 1.E10*CONVRG ERPASS = CONVRG ERROLD = ERROR XOLD = XNU YOLD = YNU C RESET FOR NEXT ITERATION 60 P1 = P2 V1 = V2 P2 = PANS V2 = VNU 70 CONTINUE C C CHECK TO SEE IF WE HAVE COME WITHIN EPSGRD OF C TARGET, IF SO ACCEPT THAT POINT IF (ERPASS.GE.EPSGRD) GO TO 80 ERROR = ERPASS XNU = XPASS YNU = YPASS GO TO 100 C SECANT METHOD DID NOT CONVERGE NORMALLY 80 OUTCOM = FAIL X2 = XNU Y2 = YNU IF (LEVEL.GE.4) WRITE (MOUTPT,99996) X2, Y2, ERROR 99996 FORMAT (15X, 38HSECANT FAILED TO CONVERGE, FOUND POINT, 2F14.9, * 14H, WITH ERROR =, 1PE11.2) RETURN C C NORMAL CONVERGENCE EXIT 90 CONTINUE C C TAKE THE BEST OF THE LAST TWO ITERATES AS THE SECANT ANSWER IF (ERROR.LE.ERROLD) GO TO 100 XNU = XOLD YNU = YOLD ERROR = ERROLD 100 CONTINUE X2 = XNU Y2 = YNU IF (LEVEL.GE.4) WRITE (MOUTPT,99995) X2, Y2, KNTSEC, PANS, ERROR 99995 FORMAT (15X, 27HSECANT CONVERGED TO (X,Y) =, 2F14.9, 4H IN , I3, * 26H ITERATIONS. PARAM,ERROR =, F11.7, 1PE11.2) RETURN END SUBROUTINE TABLGT(MOUTPT, IPAGE, NGRIDX, NGRIDY, GTYPE) TAB 10 C C *** THIS PROGRAM PRINTS THE ARRAY GTYPE GTYPE(NGRIDX,NGRIDY), C ON OUTPUT UNIT MOUTPT. IF IPAGE = 1 THEN A NEW PAGE IS C STARTED. C INTEGER GTYPE(NGRIDX,NGRIDY) C IF (IPAGE.EQ.1) WRITE (MOUTPT,99999) 99999 FORMAT (1H1, 40(1H-)/1H , 36H DOMAIN PROCESSOR OUTPUT OF GTYPE AR, * 4HRAY /1H , 40(1H-)) WRITE (MOUTPT,99998) NGRIDX, NGRIDY 99998 FORMAT (1H //1H , 4X, 58(1H+)/1H , 4X, 1H+, 56X, 1H+/1H , 4X, 1H+, * 4X, 32HTABLE OF THE POINT TYPES ON, I4, 2H X, I4, 6H GRID, * 4X, 1H+/1H , 4X, 1H+, 56X, 1H+/1H , 4X, 58(1H+)///1H , 9X, * 30HTHE POINT XGRID(1), YGRID(1) , 22H IS AT THE LOWER LEFT./1H , * 9X, 50(1H-)/) C DO 10 NY=1,NGRIDY J = NGRIDY - NY + 1 WRITE (MOUTPT,99997) J, (GTYPE(I,J),I=1,NGRIDX) 99997 FORMAT (1H , I3, 3H * , 18I7/(1H , 6X, 18I7)) 10 CONTINUE C WRITE (MOUTPT,99996) (I,I=1,NGRIDX) 99996 FORMAT (1H , 4X, 1H*/1H , 4X, 7H*******, I2, 17I7/(1H , 6X, 18I7)) C RETURN END C PROGRAM DOMTST C ******** THIS IS A MAIN PROGRAM TO TEST THE DOMAIN PROCESSOR C IT WORKS WITH PROGRAMS DOMN20 AND SETRNG WHICH DEFINE A SET C OF TWENTY TEST DOMAINS - MOST ARE PARAMETERIZED C DOMN20 - GIVES THE BOUNDARY IN PARAMETERIZED FORM C X = X(P) Y = Y(P) C SETRNG - SETS THE RANGES OF THE PARAMETERS FOR THE PIECES C AND INITIALIZES OTHER VARIABLES C THE 20 DOMAINS HAVE DEFAULF PARAMETER VALUES IN THE ARRAYS C ADFALT(20) BDFALT(20) C WHICH ARE USED IF THE SWITCH DEFALT IS .TRUE. C THE CASES ARE SELECTED BY FREE FORMAT READ DEFINED BELOW C PARAMETER (NGDIMX=43) PARAMETER (NGDIMY=43) PARAMETER (NBDIM=400) PARAMETER (NPDIM=20) C DECLARATIONS FOR 20 TEST CASES EXTERNAL DOMN20 REAL ADFALT(20),BDFALT(20) LOGICAL DEFALT INTEGER DOMAIN C COMMON / SELECT / DOMAIN,DEFALT,ADFALT,BDFALT,APAR,BPAR,MOUTPT C DECLARATIONS FOR DOMAIN PROCESSOR - USES FORTRAN 77 REAL XGRID(NGDIMX),YGRID(NGDIMY),BRANGE(2,NPDIM) DIMENSION GTYPE(NGDIMX, NGDIMY),XBOUND(NBDIM),YBOUND(NBDIM), A BPTYPE(NBDIM),BNEIGH(NBDIM),BGRID(NBDIM),BPARAM(NBDIM), B PIECE(NBDIM) INTEGER GTYPE,PIECE,BGRID,BPTYPE,BNEIGH LOGICAL FATAL,CLOCKW C THE DEFAULT PARAMETER VALUES DATA ADFALT / A .8,1.,2.,1.,.7,0.,.4,1.,.2,1.,1.,4.4,.25,0.,.5 ,1.25,3*0.,.5/ DATA BDFALT / A .6,1.,1.,1.,.1,0.,.2,1.,.2,1.,0.,3.61,1.,0.,1.5, .75,3*0.,1./ CASE = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17-19 20 C THE OUTPUT UNIT IS ALSO SET IN REGION TO 6 MOUTPT = 6 MINPUT = 5 C READ DEFAULT SWITCH READ(MINPUT,5) DEFALT 5 FORMAT(L2) C FREE FORMAT INPUT OF DATA FOR TEST CASES C INPUT VARIABLES ARE C CASE = NUMBER OF TEST DOMAIN C XA,XB = X-COORDINATES OF GRID C NXM1 = NUMBER OF X GRID SQUARES C YA,YB = Y-COORDINATES OF GRID C NYM1 = NUMBER OF Y GRID SQUARES C LEVEL = OUTPUT LEVEL FOR DOMAIN PROCESSOR C APAR, BPAR = PARAMETERS OF TEST DOMAIN C SEE SETRNG FOR SOME RESTRICTIONS ON PARAMETERS C 10 READ(MINPUT,*,END=40) A DOMAIN,XA,XB,NXM1,YA,YB,NYM1,LEVEL,APAR,BPAR NGRIDX = NXM1+1 NGRIDY = NYM1+1 C CHECK GRID SIZE IF( NGRIDX .LE. NGDIMX ) GO TO 18 WRITE(MOUTPT,12) NGDIMX 12 FORMAT(///10X,5(3H***),19H FATAL INPUT ERROR / A 10X,38HTEST DRIVER DIMENSION LIMITS X-GRID TO,I4,6H LINES) GO TO 10 18 CONTINUE IF( NGRIDY .LE. NGDIMY ) GO TO 19 WRITE(MOUTPT,13) NGDIMY 13 FORMAT(///10X,5(3H***),19H FATAL INPUT ERROR / A 10X,38HTEST DRIVER DIMENSION LIMITS Y-GRID TO,I4,6H LINES) GO TO 10 19 CONTINUE DX = (XB-XA)/NXM1 DY = (YB-YA)/NYM1 DO 20 IX = 1,NGRIDX XGRID(IX) = XA+(FLOAT(IX)-1.0)*DX 20 CONTINUE DO 30 IY = 1,NGRIDY YGRID(IY) = YA+(FLOAT(IY)-1.0)*DY 30 CONTINUE C CALL SETRNG(NBOUND,BRANGE,CLOCKW) C WRITE(MOUTPT,15) XA,XB,NXM1,YA,YB,NYM1,LEVEL 15 FORMAT(//9X,7HX-GRID ,2F10.5,6H WITH ,I2,10H INTERVALS A /9X,7HY-GRID ,2F10.5,6H WITH ,I2,10H INTERVALS B /9X,15HOUTPUT LEVEL = ,I2 ) C CALL REGION(XGRID,YGRID,NGRIDX,NGRIDY,BRANGE,NBOUND,DOMN20, A CLOCKW,.FALSE.,LEVEL, B GTYPE,XBOUND,YBOUND,PIECE,BPTYPE,BNEIGH, C BGRID,BPARAM,NBDIM,NBPTS,FATAL) IF( FATAL ) WRITE(MOUTPT,35) 35 FORMAT(42H ++ ++ ++ FATAL ERROR IN DOMAIN PROCESSING) GO TO 10 C 40 CONTINUE WRITE(MOUTPT,45) 45 FORMAT(//,20X,24HEND OF DOMAIN PROCESSING) STOP END SUBROUTINE DOMN20(P,X,Y,IPIECE) C C ***** THIS IS A USER SUPPLIED ROUTINE TO DEFINE THE DOMAIN BOUNDARY C THIS ONE HAS 20 CASES THAT CAN BE SELECTED C MOST DOMAINS HAVE 2 PARAMETERS A AND B TO VARY THEM. C DEFAULT VALUES ARE IN THE ADFALT, BDFALT ARRAYS SET IN TEST20 C THE DEFAULT VALUES ARE LISTED HERE FOR REFERENCE C DATA ADFALT / C A .8,1.,2.,1.,.7,0.,.4,1.,.2,1.,1.,4.4,.25,0.,.5 ,1.25,3*0.,.5/ C DATA BDFALT / C A .6,1.,1.,1.,.1,0.,.2,1.,.2,1.,0.,3.6, 1.,0.,1.5, .75,3*0.,1./ CASE = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17-19 20 C REAL ADFALT(20),BDFALT(20) LOGICAL DEFALT INTEGER DOMAIN COMMON / SELECT / DOMAIN,DEFALT,ADFALT,BDFALT,APAR,BPAR,MOUTPT DATA PI A / 3.14159265359 / C IF( DEFALT ) THEN A = ADFALT(DOMAIN) B = BDFALT(DOMAIN) ELSE A = APAR B = BPAR ENDIF GO TO (100,200,300,400,500,600,700,800,900,1000,1100,1200,1300, A 1400,1500,1600,1700,1800,1900,2000) ,DOMAIN C C DOMAIN 1 = RECTANGLE WITH CORNER 0,0 C PARAMTERS = COORDINATES OF UPPER RIGHT CORNER C NBOUND = 4 C BRANGE 1 = 0.0 0.0 0.0 0.0 C BRANGE 2 = B A B A 100 GO TO( 110,120,130,140),IPIECE 110 Y= P X= 0.0 RETURN 120 X = P Y = B RETURN 130 Y= B - P X= A RETURN 140 X = A - P Y = 0.0 RETURN C C DOMAIN 2 = TRIANGLE WITH BASE 0 TO 1 C PARAMETERS = THIRD VERTEX COORDINATES C NBOUND = 3 C BRANGE 1 = 0.0 0.0 0.0 C BRANGE 2 = B B 1.0 200 CONTINUE GO TO( 210,220,230),IPIECE 210 Y = P X = A*Y/B RETURN 220 Y = B - P X = (A-1.)*Y/B + 1. RETURN 230 X = 1.0 - P Y = 0.0 RETURN C C DOMAIN 3 = ELLIPSE, COUNTERCLOCKWISE C PARAMETERS = MAJOR, MINOR SEMI-AXES C NBOUND = 1 BRANGE 1 = 0.0 BRANGE 2 = 2.*PI 300 CONTINUE X = A*COS(P) Y = B*SIN(P) RETURN C C DOMAIN 4 = SEVEN-SIDED OBJECT C PARAMETERIZED BY THE MAPPING C X = U**A Y = V**B C NBOUND = 7 C BRANGE 1 = 0.0 0.2 0.6 1.2 0.0 0.0 0.0 C BRANGE 2 = 0.2 0.6 1.2 1.4 0.4 0.5 0.5 C 400 GO TO(410,420,430,440,450,460,470),IPIECE 410 U = P V = 5.*U GO TO 499 420 U = P V = 1.2 - U GO TO 499 430 U = P V = 0.666666666667*U + .2 GO TO 499 440 U = P V = -5.*U + 7.0 GO TO 499 450 U = 1.4 - P V = 0.5*P GO TO 499 460 U = 1.0 - P V = 0.2 - .4*P GO TO 499 470 U = 0.5 - P V = 0.0 499 X = U**A Y = V**B RETURN C C DOMAIN 5 = EIGHT SIDED OBJECT C PARAMETERS = COORDINATES OF LAST VERTEX 500 CONTINUE GO TO(510,520,530,540,550,560,570,580),IPIECE 510 X = .08 + P Y = .7 + 1.5625*P RETURN 520 X = .4 + P Y =1.2 - 1.3666666666667*P RETURN 530 X = .7 + P Y = .79 + 4.*P RETURN 540 X = .85 + P Y =1.39 - 1.2*P RETURN 550 X =1.35 + P Y = .79 + .6*P RETURN 560 X =1.8 Y =1.06 - P RETURN 570 S = (.3-B)/(1.8-A) X =1.8 - P Y = S*(X-1.8) + .3 RETURN 580 S = (.7-B)/(.08-A) X = A - P Y = S*(X-A) + B RETURN C C DOMAIN 6 = SIXTEEN SIDES WITH OSCILATIONS ABOUT THE C LINES Y = 5 AND Y = 1 C PARAMETERIZED BY THE MAPPING C X = U + A*LOG(1+V) + B*SQRT(U*V) C Y = V - B/SQRT(1+U+V) + A*EXP((U-V)/(U+V)) 600 CONTINUE GO TO(601,602,603,604,605,606,607,608,609,610, A 611,612,613,614,615,616), IPIECE 601 U = 9.2 - .1481481481481*P V = 5.2 - P*17./27. GO TO 699 602 U = 8.8 + .2592592592592*P V = 3.5 - .6666666666667*P GO TO 699 603 U = 9.5 - P V = 1.7 - .833333333333*P GO TO 699 604 U = 8.3 - P V = .7 + P GO TO 699 605 U = 7. - P*2.2 V = 2. - P*1.4 GO TO 699 606 U = 4.8 - P V = 0.6 + P GO TO 699 607 U = 4.0 - P V = 1.4 - .5*P GO TO 699 608 U = 2.8 - P V = 0.8 + .222222222222*P GO TO 699 609 U = 1.45- P V = 1.1 - .2*P GO TO 699 610 U = 0.45 + P V = 0.9 + 2.2857142857*P GO TO 699 611 U = 2.2 + P V = 4.9 + P/3. GO TO 699 612 U = 3.4 + P V = 5.3 - .6*P GO TO 699 613 U = 4.4 + P V = 4.7 + 2.0*P GO TO 699 614 U = 5.5 + P V = 6.9 - 1.166666666667*P GO TO 699 615 U = 7.3 + P V = 4.8 + .57142857143*P GO TO 699 616 U = 8.0 + P V = 5.2 699 X = U + A*ALOG(1.+V) + B*SQRT(U*V) Y = V - B/SQRT(1.+U+V) + A*EXP((U-V)/(U+V)) RETURN C C DOMAIN 7 = HALF OF SQUARE WITH CIRCULAR HOLE, NBOUND = 6 C PARAMETERS = LOCATION OF CENTER ON X-AXIS AND RADIUS 700 GO TO(710,720,730,740,750,760), IPIECE 710 X = P Y = .4 RETURN 720 X = .8 Y = .4 - P RETURN 730 X = .8 - P Y = .0 RETURN 740 X = A + B*COS(PI*P) Y = B*SIN(PI*P) RETURN 750 X = A-B - P Y = .0 RETURN 760 X = .0 Y = P RETURN C C DOMAIN 8 = CIRCULAR DISK C PARAMETERS ARE CENTER (A,A) AND RADIUS B 800 CONTINUE X = A - B*COS(PI*P) Y = A + B*SIN(PI*P) RETURN C C DOMAIN 9 = MEAT CLEAVER WITH WAVY TOP, NBOUND = 6 C ********* THIS DOMAIN REQUIRES A VERY FINE GRID ************ C PARAMETERS = AMPLITUDES OF 2 WAVES AT TOP LEFT 900 GO TO(910,920,930,940,950,960), IPIECE 910 X = P IF( P .GE. 3 ) GO TO 911 Y = 5. - A*SIN(PI*P) RETURN 911 IF( P .GE. 4. ) GO TO 912 Y = 5. + B*SIN(2.*PI*(P-3.)) RETURN 912 IF( P .GE. 5.5) GO TO 913 Y = 5. + .15*SIN(1.33333333333*PI*(P-4.)) RETURN 913 IF( P .GE. 6.2) GO TO 914 Y = 5. + .15*SIN(2.85714285714*PI*(P-5.5)) RETURN 914 Y = 5. + .12*SIN(5.*PI*(P-6.2)) RETURN 920 X = 7.5 Y = 5. + .12*SIN(6.5*PI) - P RETURN 930 X = 7.5 - P P930 = AMIN1(AMAX1(P,-2.0),6.0) Y = 1.1 + EXP(-2.*P930) RETURN 940 X = 3.5 - P Y = 1.1003354626279 + 1.900896383536*P**3*(8.+(P-.5)**3) RETURN 950 X = 3.0 - P Y = 3.001231846164 RETURN 960 X = - .99938407691*SIN(PI*P) Y = 3.001231846164 + 1.998768153836*P RETURN C C DOMAIN 10 = TRIANGLE WITH ONE CORNER SMOOTHED, NBOUND = 4 C PARAMETERS = RADIUS OF SMOOTHING ( ONLY 1 PARAMETER ) 1000 CONTINUE GO TO(1010,1020,1030), IPIECE 1010 X = -P Y = .75*P RETURN 1020 X =-4.0 + P Y = 3.0 RETURN 1030 IF( P .GT. 0.5 ) GO TO 1040 X = -A + A*SIN(PI*P) Y = 3.0-A + A*COS(PI*P) RETURN 1040 X = 0.0 Y = 3.5 -A - P RETURN C C DOMAIN 11 = ALMOST BOX WITH CURVED RIGHT SIDE, NBOUND = 4 C PARAMETERIZED BY MAPPING DOMAIN WITH C X = A*U/(1+B*V) Y = COS(B)*V/(1+.5*SIN((1-A)*U)) 1100 GO TO(1110,1120,1130,1140), IPIECE 1110 U = 1.0 V = 0.5 + P GO TO 1199 1120 U = 1.0 + P V = 4.0 GO TO 1199 1130 U = 4.0 + .1*P*(P-4.5)**2 V = 4.0 - P GO TO 1199 1140 U = 4.0 - P V=-0.5+P*.333333333333 1199 X = A*U/(1.+B*V) Y = COS(B)*V/(1.+.5*SIN((1.-A)*U)) RETURN C C DOMAIN 12 = COMPLICATED WITH LONG, NARROW NECKS. NBOUND = 14 C PARAMETERS ADJUST WIDTH AND LENGTH OF ONE LONG NECK 1200 GO TO(1201,1202,1203,1204,1205,1206,1207,1208,1209 A ,1210,1211,1212,1213,1214),IPIECE 1201 X = 1.0 + P Y = 5.4 RETURN 1202 X = 11.5 Y = 5.4 - P RETURN 1203 X = 11.5 - P Y = 4.7 RETURN 1204 X = A Y = 4.7 - P RETURN 1205 X = A + P Y = B RETURN 1206 P1206 = AMIN1(AMAX1(P,-1.),10.) X = 13.5 - COS(PI*P1206) P15 = AMIN1(P1206,1.5) Y = B + (B-3.4)*SIN(PI*P15)*(1.+16.*P15*(P15-1.5)**2) RETURN 1207 X = 12.5 - P Y = 3.4 RETURN 1208 X = 7.2 Y = 3.4 - P RETURN 1209 X = 7.2 - P Y = 0.3 RETURN 1210 P1210 = AMIN1(1.0,AMAX1(P,0.0)) X = 3.6 - 6.*SQRT(P1210)*(1.-P) Y = 0.3 + .4* P - 0.2*P*(P-2.)*SQRT(1.0-P1210) RETURN 1211 X = 3.6 + P Y = 0.7 RETURN 1212 X = 6.4 Y = 0.7 + P RETURN 1213 X = 6.4 - P Y = 2.7 RETURN 1214 X = 0.4 + P Y = 2.7 + 4.5*P RETURN C C DOMAIN 13 = BOX WITH WAVY TOP, NBOUND = 4 C PARAMETERS = AMPLITUDE AND FREQUENCY OF TOP 1300 GO TO(1310,1320,1330,1340), IPIECE 1310 X = 2.*P IB = B Y = 1.-A + A*COS(PI*P*2.*FLOAT(IB)) RETURN 1320 X = 4.0 Y = 1.0 - P RETURN 1330 X = 4.0 - P Y = 0.0 RETURN 1340 X = 0.0 Y = P RETURN C C DOMAIN 14 = PROBLEM 17 OF OLD PDE EVAL. WITH DEFAULT PARAM. C REGION FROM REACTOR WALL HEAT TRANSFER PROBLEM, NBOUND = 9 C PARAMETERS MAP DOMAIN BY C X = U + A*V*(U-.5) C Y = V + B*(1.-(1.5*U-.7)**2) 1400 GOTO(1410,1415,1420,1425,1430,1440,1450,1460,1470,1480,1490) A ,IPIECE 1410 U = P V = 1.0 GO TO 1499 1415 U = 0.65 V = 1.0 - P GO TO 1499 1420 U = 0.65 + P V = 0.7 GO TO 1499 1425 U = 0.85 V = 0.7 - P GO TO 1499 1430 U = 0.85 + P V = 0.5 GO TO 1499 1440 U = .95 V = 0.5 - P GO TO 1499 1450 U = .95 + P V = .3 GO TO 1499 1460 U = 1.0 V = .3 - P GO TO 1499 1470 U = 1.0 - P V = 0.0 GO TO 1499 1480 U = .8*COS(P*PI) V = .8*SIN(P*PI) GO TO 1499 1490 U = 0.0 V = .8 + P 1499 X = U + A*V*(U-.5) Y = V + B*(1.-(1.5*U-.7)**2) RETURN C C DOMAIN 15 = QUARTER ANNULUS, NBOUND = 4 C PARAMETERS = INNER AND OUTER RADII 1500 GO TO(1510,1520,1530,1540),IPIECE 1510 X = A*COS(PI*P) Y = A*SIN(PI*P) RETURN 1520 X = 0. Y = A + P RETURN 1530 X = B*SIN(PI*P) Y = B*COS(PI*P) RETURN 1540 X = B - P Y = 0.0 RETURN C C DOMAIN 16 = REGION FROM HALL,LUCZAK,SERDY - TOMS 2(1976)257-274 C NBOUND = 5, NEAR TRIANGLE WITH SEMI-CIRCULAR NOTCH C PARAMETERS MODIFY NOTCH 1600 GO TO(1610,1620,1630,1640,1650),IPIECE 1610 X = -2. + 4.5*COS(.4605539917-P) Y = -2. + 4.5*SIN(.4605539917-P) RETURN 1620 X = 2.5 - P Y = -2.0 RETURN 1630 X = A + B*COS(PI*P) Y =- 2.0 + B*SIN(PI*P) RETURN 1640 X = A-B - P Y = -2.0 RETURN 1650 X = -2. + 2.015564437*P Y = P - 2.0 RETURN C C DOMAIN 17 = SECOND DOMAIN FROM HALL ET AL C SPATULA SHAPE WITH NBOUND = 6 C PARAMETERIZED BY ROTATING A RADIANS, EXPANDING BY 1+B 1700 GO TO(1710,1720,1730,1740,1750,1760),IPIECE 1710 U = 0.0 V = P GO TO 1799 1720 U = P V = 23. GO TO 1799 1730 IF( P .GE. 5. ) GO TO 1735 P1730 = AMAX1(AMIN1(P,1000.),-2.) U = 20. - 10.*COS(PI*P1730*.1) V = 23. - 10.*SIN(PI*P1730*.1) GO TO 1799 1735 U = 15. + P V = 13. GO TO 1799 1740 U = 34. V = 13. - P GO TO 1799 1750 IF( P .GT. 15. ) GO TO 1755 U = 34. - P V = 7. GO TO 1799 1755 P1755 = AMAX1(AMIN1(P,75.),-25.) U = 19. - 7.*SIN(PI*(P1755-15.)*.1) V = 7.*COS(PI*(P1755-15.)*.1) GO TO 1799 1760 U = 12. - P V = 0.0 1799 X = (1.+B)*( U*COS(A) - V*SIN(A) ) Y = (1.+B)*( V*COS(A) + U*SIN(A) ) RETURN C C DOMAIN 18 HAS LOTS OF TANGENTS TO LINES IN DEFAULT CASE C X = -1. 1.3 3.0 NBOUND = 7 C Y = 0. 0.5 -.5 C BRANGE = 0 0 0 -1 -S 0 -S C (S=.6*PI) S 1 S 1 0 (7-SQRT(3))/2 PI/2 C PARAMETERS MAP DOMAIN BY C X = U + .1*SIN(V*A) C Y = V + B*U**2 1800 GO TO(1810,1820,1830,1840,1850,1860,1870), IPIECE 1810 U = -1. - COS(P) V = SIN(P) GO TO 1899 1820 XS = -1. - COS(.6*PI) YS = SIN(.6*PI) U = XS + P*(1.5-XS) V = YS + P*(1.5-YS) GO TO 1899 1830 U = 2.0 - .5*COS(P) V = 1.5 - .5*SIN(P) GO TO 1899 1840 XS = 2. - .5*COS(.6*PI) YS = 1.5- .5*SIN(.6*PI) V = -P*YS U = 3. -(3.-XS)*P**2 GO TO 1899 1850 U = 2.0 - .5*COS(P) V =-1.5 - .5*SIN(P) GO TO 1899 1860 U = 1.5 - P V = -1.5 GO TO 1899 1870 U = -2. + COS(P) V = -1. + SIN(P) 1899 X = U + .1*SIN(V*A) Y = V + B*U**2 RETURN C C DOMAIN 19 = COMPLEX REGION WITH SHARP INTERIOR REENTRANT POINT C ORINTATION IS COUNTER CLOCKWISE, NBOUND = 7 C PARAMETERS MAP DOMAIN BY C X = U*(1.+A*U**2*(V-2.)**2/10.) C Y = V*(1.+B*U**2*(V-1.)**2/10.) C ******** THIS DOMAIN REQUIRES A VERY FINE GRID *************** 1900 GO TO(1910,1920,1930,1940,1950,1960,1970), IPIECE 1910 U = 2.*SIN(P*PI) V = 2. -2.*COS((P+P**2/4.)*PI) GO TO 1999 1920 XS = 2.*SIN(.6*PI) YS = 2. - 2.*COS(.69*PI) U = XS - .9*P V = YS - .8*P -P*(1.-P) GO TO 1999 1930 XS = 2.*SIN(.6*PI) - .9 YS = 2. - 2.*COS(.69*PI) - .8 U = P + (1.-P)*XS V = P + (1.-P)*YS GO TO 1999 1940 U = 1. - P/100. V = 1. + 1.4*P GO TO 1999 1950 XS = 2.*SIN(.64*PI) YS = 2. - 2.*COS(.7424*PI) U = 1.07 + P*(XS-.99) - .16*ABS(.5-P) V = 2.4 + P*(YS-2.4) GO TO 1999 1960 U = 2.*SIN(P*PI) V = 2.-2.*COS(PI*(P+P**2/4.*(1.64-P))) GO TO 1999 1970 U = -2.*P*(1.-P) YS = 2.-2.*COS(PI*1.64) V = YS*(1.-P)-P*(1.-P)*(1.+2.*P*(1.-P)+SIN(6.*P*PI)*2.*ABS(P-.5)) 1999 X = U*(1.+A*U**2*(V-2.)**2/10.) Y = V*(1.+B*U**2*(V-1.)**2/10.) RETURN C C DOMAIN 20 = TWO INTERSECTING CIRCLES, NBOUND = 2 C CIRCLES INTERSECT AT X=A AND ARE COUNTERCLOCKWISE C LEFT ONE HAS RADIUS 1., RIGHT ONE HAS RADIUS B 2000 GO TO(2010,2020),IPIECE 2010 X = A + SQRT(A*A+B*B-1.) + B*COS(P) Y = B*SIN(P) RETURN 2020 X = COS(P) Y = SIN(P) RETURN END SUBROUTINE SETRNG(NBOUND,BRANGE,CLOCKW) C C **** THIS PROGRAM SETS THE BOUNDARY DATA FOR THE KASE AT HAND REAL BRANGE(2,25),ADFALT(20),BDFALT(20) INTEGER DOMAIN LOGICAL CLOCKW,DEFALT COMMON / SELECT / DOMAIN,DEFALT,ADFALT,BDFALT,APAR,BPAR,MOUTPT DATA PI ,LEVEL A / 3.14159265359 , 2 / C SET PARAMETERS FOR DOMAIN IF( DEFALT ) THEN A = ADFALT(DOMAIN) B = BDFALT(DOMAIN) ELSE A = APAR B = BPAR ENDIF CLOCKW = .TRUE. C SELECT DOMAIN GO TO (100,200,300,400,500,600,700,800,900,1000,1100,1200,1300, A 1400,1500,1600,1700,1800,1900,2000) ,DOMAIN C C DOMAIN 1 = RECTANGLE( A,B .GT. 0.0 ) 100 NBOUND = 4 BRANGE(1,1) = 0.0 BRANGE(1,2) = 0.0 BRANGE(1,3) = 0.0 BRANGE(1,4) = 0.0 BRANGE(2,1) = B BRANGE(2,2) = A BRANGE(2,3) = B BRANGE(2,4) = A GO TO 2100 C DOMAIN 2 = TRIANGLE( B .GT. 0.0 ) 200 NBOUND = 3 BRANGE(1,1) = 0.0 BRANGE(1,2) = 0.0 BRANGE(1,3) = 0.0 BRANGE(2,1) = B BRANGE(2,2) = B BRANGE(2,3) = 1.0 GO TO 2100 C DOMAIN 3 = ELLIPSE 300 NBOUND = 1 CLOCKW = .FALSE. BRANGE(1,1) = 0.0 BRANGE(2,1) = 2.*PI GO TO 2100 C DOMAIN 4 = SEVEN-SIDED OBJECT 400 NBOUND = 7 BRANGE(1,1) = 0.0 BRANGE(1,5) = 0.0 BRANGE(1,6) = 0.0 BRANGE(1,7) = 0.0 BRANGE(1,2) = 0.2 BRANGE(2,1) = 0.2 BRANGE(2,6) = 0.5 BRANGE(2,7) = 0.5 BRANGE(1,3) = 0.6 BRANGE(2,2) = 0.6 BRANGE(1,4) = 1.2 BRANGE(2,3) = 1.2 BRANGE(2,5) = 0.4 BRANGE(2,4) = 1.4 GO TO 2100 C DOMAIN 5 = EIGHT SIDED OBJECT( .08 LT A LT 1.8 ) 500 NBOUND = 8 DO 510 J=1,8 510 BRANGE(1,J) = 0.0 BRANGE(2,1) = .32 BRANGE(2,2) = .3 BRANGE(2,3) = .15 BRANGE(2,4) = .5 BRANGE(2,5) = .45 BRANGE(2,6) = .76 BRANGE(2,7) = 1.8 - A BRANGE(2,8) = A - .08 GO TO 2100 C DOMAIN 6 = SIXTEEN SIDES WITH OSCILATIONS ABOUT THE C LINES X = 5 AND X = 1 600 NBOUND = 16 DO 610 J = 1,16 610 BRANGE(1,J) = 0.0 BRANGE(2,1 )= 2.7 BRANGE(2,2 )= 2.7 BRANGE(2,3 )= 1.2 BRANGE(2,4 )= 1.3 BRANGE(2,5 )= 1.0 BRANGE(2,6 )= 0.8 BRANGE(2,7 )= 1.2 BRANGE(2,8 )= 1.35 BRANGE(2,9 )= 1.0 BRANGE(2,10)= 1.75 BRANGE(2,11)= 1.2 BRANGE(2,12)= 1.0 BRANGE(2,13)= 1.1 BRANGE(2,14)= 1.8 BRANGE(2,15)= 0.7 BRANGE(2,16)= 1.2 GO TO 2100 C DOMAIN 7 = HALF OF SQUARE WITH CIRCULAR HOLE C (A+B) .LT. .8 A .GT. B 700 NBOUND = 6 DO 710 J = 1,6 710 BRANGE(1,J) = 0.0 BRANGE(2,1) = 0.8 BRANGE(2,2) = 0.4 BRANGE(2,3) = .8 - (A+B) BRANGE(2,4) = 1.0 BRANGE(2,5) = A - B BRANGE(2,6) = 0.4 GO TO 2100 C DOMAIN 8 = CIRCULAR DISK 800 NBOUND = 1 BRANGE(1,1) = 0.0 BRANGE(2,1) = 2.0 GO TO 2100 C DOMAIN 9 = MEAT CLEAVER WITH WAVY TOP 900 NBOUND = 6 DO 910 J = 1,6 910 BRANGE(1,J) = 0.0 BRANGE(2,1) = 7.5 BRANGE(2,2) = 2.9 + .12*SIN(6.5*3.14159265359) BRANGE(2,3) = 4.0 BRANGE(2,4) = 0.5 BRANGE(2,5) = 3.0 BRANGE(2,6) = 1.0 GO TO 2100 C DOMAIN 10 = TRIANGLE WITH ROUNDED CORNER C 0.0 .LT. A .LT. 3. 1000 NBOUND = 3 BRANGE(1,1) = 0.0 BRANGE(1,2) = 0.0 BRANGE(1,3) = 0.0 BRANGE(1,4) = 0.0 BRANGE(2,1) = 4.0 BRANGE(2,2) = ( 4. - A ) BRANGE(2,3) = 3.5 - A GO TO 2100 C DOMAIN 11 = ALMOST BOX, CURVED RIGHT SIDE 1100 NBOUND = 4 BRANGE(1,1) = 0.0 BRANGE(1,2) = 0.0 BRANGE(1,3) = 0.0 BRANGE(1,4) = 0.0 BRANGE(2,1) = 3.5 BRANGE(2,2) = 3.0 BRANGE(2,3) = 4.5 BRANGE(2,4)=3.0 GO TO 2100 C DOMAIN 12 = COMPLICATED WITH LONG, NARROW NECKS 1200 NBOUND = 14 DO 1210 J = 1,14 1210 BRANGE(1,J) = 0.0 BRANGE(2,1) = 10.5 BRANGE(2,2) = 0.7 BRANGE(2,3) = 11.5 - A BRANGE(2,4) = 4.7 - B BRANGE(2,5) = 12.5 - A BRANGE(2,6) = 2.0 BRANGE(2,7) = 5.3 BRANGE(2,8) = 3.1 BRANGE(2,9) = 3.6 BRANGE(2,10)= 1.0 BRANGE(2,11)= 2.8 BRANGE(2,12)= 2.0 BRANGE(2,13)= 6.0 BRANGE(2,14)= 0.6 GO TO 2100 C DOMAIN 13 = BOX WITH WAVY TOP 1300 NBOUND = 4 BRANGE(1,1) = 0.0 BRANGE(1,2) = 0.0 BRANGE(1,3) = 0.0 BRANGE(1,4) = 0.0 BRANGE(2,1) = 2.0 BRANGE(2,2) = 1.0 BRANGE(2,3) = 4.0 BRANGE(2,4) = 1.0 GO TO 2100 C C DOMAIN 14 = PROBLEM 17 FROM PDE EVALUATION 1400 NBOUND = 11 DO 1410 J = 1,NBOUND 1410 BRANGE(1,J) = 0.0 BRANGE(2,1) = .65 BRANGE(2,2) = .3 BRANGE(2,3) = .2 BRANGE(2,4) = .2 BRANGE(2,5) = .1 BRANGE(2,6) = .2 BRANGE(2,7) = .05 BRANGE(2,8) = .3 BRANGE(2,9) = .2 BRANGE(2,10) = .5 BRANGE(2,11) = .2 GO TO 2100 C DOMAIN 15 = QUARTER ANNULUS( A .LT. B ) 1500 NBOUND = 4 DO 1510 J = 1,NBOUND 1510 BRANGE(1,J) = 0.0 BRANGE(2,1) = 0.5 BRANGE(2,2) = B - A BRANGE(2,3) = 0.5 BRANGE(2,4) = B - A GO TO 2100 C DOMAIN 16 = ALMOST TRIANGLE WITH NOTCH C (A+B) .LT. 2.5 B .LT. A 1600 NBOUND = 5 DO 1610 J = 1,NBOUND 1610 BRANGE(1,J) = 0.0 BRANGE(2,1) = 0.4605539917 BRANGE(2,2) = 2.5 - (A + B) BRANGE(2,3) = 1.0 BRANGE(2,4) = 2.0 + A-B BRANGE(2,5) = 2.0 GO TO 2100 C DOMAIN 17 = SPATULA SHAPE 1700 NBOUND = 6 DO 1710 J = 1,NBOUND 1710 BRANGE(1,J) = 0.0 BRANGE(2,1) = 23.0 BRANGE(2,2) = 10.0 BRANGE(2,3) = 19.0 BRANGE(2,4) = 6.0 BRANGE(2,5) = 20.0 BRANGE(2,6) = 12.0 GO TO 2100 C DOMAIN 18 HAS LOTS OF TANGENTS 1800 NBOUND = 7 DO 1810 I = 1,7 1810 BRANGE(1,I) = 0 YS = .6*PI BRANGE(1,4) = -1.0 BRANGE(1,5) = -YS BRANGE(1,7) = -PI/6. BRANGE(2,1) = YS BRANGE(2,2) = 1.0 BRANGE(2,3) = YS BRANGE(2,4) = 1.0 BRANGE(2,5) = 0.0 BRANGE(2,6) = (7.-SQRT(3.))/2. BRANGE(2,7) = PI/2. GO TO 2100 C DOMAIN 19 = VERY COMPLICATED DOMAIN 1900 NBOUND = 7 CLOCKW = .FALSE. DO 1910 I = 1,7 1910 BRANGE(1,I) = 0.0 BRANGE(1,6) = 0.64 BRANGE(2,1) = 0.6 BRANGE(2,2) = 1.0 BRANGE(2,3) = 1.0 BRANGE(2,4) = 1.0 BRANGE(2,5) = 1.0 BRANGE(2,6) = 2.0 BRANGE(2,7) = 1.0 GO TO 2100 C DOMAIN 20 = TWO INTERSECTING DISKS C A .LT. 1. SQRT(1-A*A) .LT. B 2000 CONTINUE NBOUND = 2 CLOCKW = .FALSE. ANGLE1 = ATAN(1./SQRT(B*B/(1.-A*A)-1.)) ANGLE2 = ATAN(SQRT(1.-A*A)/A) BRANGE(1,1) = ANGLE1 - PI BRANGE(2,1) = -BRANGE(1,1) BRANGE(1,2) = ANGLE2 BRANGE(2,2) = 2.*PI - ANGLE2 C 2100 IF( LEVEL .LT. 2 ) GO TO 2150 C PRINT PROBLEM PARAMETERS WRITE(MOUTPT,2125) DOMAIN,NBOUND, A (BRANGE(1,I),BRANGE(2,I),I=1,NBOUND) 2125 FORMAT(1H1,/,9X,13HDOMAIN NUMBER,I3, A 06H WITH ,I2,16H BOUNDARY PIECES,/, A 21X,21HBRANGE PARAMETERS ARE,/,(15X,2F12.7)) IF( .NOT. CLOCKW ) WRITE(MOUTPT,2127) 2127 FORMAT(9X,29HBOUNDARY IS COUNTER-CLOCKWISE) IF( DEFALT ) THEN WRITE(MOUTPT,2130) A,B 2130 FORMAT(9X,29HDEFAULT DOMAIN PARAMETERS ARE, 2F12.6) ELSE WRITE(MOUTPT,2135) A,B 2135 FORMAT(9X,21HDOMAIN PARAMETERS ARE,2F12.6) END IF 2150 CONTINUE RETURN END T THIS IS TEST DATA SET 1. 1 0.0 0.8 3 0.0 1.0 5 2 0. 0. 1 0.0 1.0 5 0.0 1.0 5 2 0. 0. 2 0.0 1.0 5 0.0 1.0 5 2 0. 0. 3 -2.0 2.0 5 -2.0 2.0 5 2 0. 0. 4 0.0 1.5 10 -0.1 1.1 8 2 0. 0. 5 0.0 2.0 12 0.0 1.5 10 2 0. 0. 6 0.0 10.0 20 -1.0 7.0 16 2 0. 0. 7 0.0 1.0 20 0.0 1.0 10 2 0. 0. 8 -0.5 3.0 6 -0.5 3.0 6 2 0. 0. 9 -2.0 8.0 10 1.0 6.0 20 2 0. 0. 10 -5.0 1.0 6 0.0 3.0 3 2 0. 0. 11 -1.0 6.0 7 -1.0 6.0 7 2 0. 0. 12 0.0 15.0 30 0.0 6.0 20 2 0. 0. 13 0.0 5.0 10 -1.1 1.1 5 2 0. 0. 14 0.0 1.0 25 0.0 1.05 14 2 0. 0. 15 0.0 2.0 10 0.0 2.0 10 2 0. 0. 16 -2.5 2.5 10 -2.0 0.0 8 2 0. 0. 17 0.0 36.0 12 0.0 24.0 8 2 0. 0. 18 -2.0 4.0 6 -1.5 2.0 7 2 0. 0. 19 -2.0 2.0 15 -0.4 4.0 11 2 0. 0. 20 -1. 2.6 12 -1.1 1.1 11 2 0. 0. F THIS IS TEST DATA SET 2. 1 0.0 1.0 5 0.0 1.0 5 2 .9 .3 2 0.0 1.0 5 0.0 1.0 5 2 .5 .8 3 -2.0 2.0 10 -2.0 2.0 7 2 1. 2. 4 0.0 1.5 10 0.0 1.0 10 2 .8 1.1 5 0.0 2.0 10 0.0 1.5 10 2 1.2 .4 6 0.0 20.0 20 -1.0 7.5 17 2 .2 .1 7 0.0 1.0 20 0.0 1.0 10 2 .3 .1 8 -0.5 3.0 7 -0.5 3.0 7 2 .8 1.1 9 -2.0 8.0 10 1.0 6.0 20 2 0.1 0. 10 -5.0 1.0 6 0.0 3.0 3 2 1.5 0.0 11 -1.0 6.0 7 -1.0 6.0 7 2 .8 .2 12 0.0 15.0 30 0.0 10.5 21 2 10. 4.1 13 0.0 5.0 20 -0.4 1.1 15 2 0.2 2. 14 -.25 1.25 25 0.0 1.4 14 2 .1 .2 15 0.0 2.0 10 0.0 2.0 10 2 .4 .9 16 -2.5 2.5 10 -2.0 0.0 8 2 1. .6 17 -20. 16.0 12 0.0 32.0 8 2 1. -.2 18 -2.5 3.5 20 -2.0 6.0 24 2 0.0 0.5 19 -4.0 4.0 20 -0.4 4.0 11 2 0.5 0.0 20 -1.0 2.5 10 -1. 1.0 11 2 .4 .92 F THIS IS TEST DATA SET 3. 1 0.0 1.0 10 0.0 .32 8 2 1.0 .1 2 0.0 2.0 10 0.0 .49 7 2 2.0 .45 3 -2.0 2.0 20 -2.0 2.0 10 2 2. .2 4 0.0 1.5 10 0.0 1.0 10 2 .8 1.1 5 0.0 2.0 10 0.0 1.5 15 2 1.7 .5 6 1.0 18.0 34 0.0 8.0 32 2 0.5 1. 7 0.0 1.0 20 0.0 .48 12 2 .3 .25 8 -0.5 1.0 15 0.0 1.0 40 2 .6 .05 9 -2.0 8.0 20 1.0 7.3 42 2 1. 1.4 10 -5.0 1.0 6 0.0 3.0 3 2 2.7 0.0 11 0.0 10.0 20 -1.0 8.0 18 2 2. .3 12 0.0 15.0 30 0.0 8.4 21 2 2.5 3.7 13 0.0 5.0 20 0.0 1.5 15 2 0.35 3. 14 -.5 1.5 25 0.0 2.0 20 2 1.0 1.0 15 0.0 2.0 20 0.0 2.0 20 2 .36 .51 16 -2.3 2.5 12 -2.0 0.0 20 2 .2 .75 17 0.0 60.0 12 -50. 30.0 8 2 -1. 1.0 18 -4.0 4. 16 -2. 4. 6 2 3. 0.0 19 -3.0 3. 15 -0.8 6.4 18 2 .2 .2 20 -1.0 3.9 14 -1.6 1.6 16 2 .8 1.5 F THIS IS TEST DATA SET 4. 1 0.0 1.0 5 0.0 .20 25 2 1.0 .01 2 -2. 1.0 9 -0.01 0.24 10 2 -2. .04 3 -0.2 0.2 10 -2.0 2.0 5 2 .02 2.0 4 0.0 3.0 30 0.0 1.0 20 2 3. 0.4 5 0.0 2.0 10 0.0 1.5 10 2 .85 .3 6 0.0 22.0 22 -3.0 7.6 16 2 -1. 1. 7 0.0 1.0 20 -0.01 0.99 25 2 .4 .36 8 -1.0 2.0 15 -0.8 1.8 13 2 .222 .888 9 -2.0 8.0 25 1.0 8.0 28 2 -1.6 1.5 10 -5.0 1.0 6 0.0 3.2 40 2 0.1 0.0 11 -1.0 20.0 21 -2.0 8.0 20 2 3. 0.5 12 0.0 15.0 30 -0.05 6.45 26 2 1.2 3.5 13 0.0 4.4 22 0.0 1.2 12 2 0.4 1.7 14 -1. 2. 20 0.0 1.5 15 2 2.0 0.5 15 0.0 1.0 20 0.0 1.0 20 2 .72 .8 16 -2.5 2.5 10 -2.0 0.0 8 2 1.5 .9 17 -32. 4.0 12 0.0 30.0 10 2 1.56 -0.25 18 -4. 4. 32 -2. 3. 25 2 1. 0.2 19 -3 3. 15 -0.8 7.2 16 2 .1 .5 20 -1.0 3.0 16 -1. 1.0 10 2 .9 1.0 C PROGRAM DOMTS2 EXTERNAL CIRCLE, CIR2,CIR3 PARAMETER( NGRID = 32 ) C REAL XGRID(NGRID),YGRID(NGRID),BRANGE(2,5) DIMENSION GTYPE(NGRID, NGRID),XBOUND(300),YBOUND(300), A BPTYPE(300),BNEIGH(300),BGRID(300),BPARAM(300), B PIECE(300) INTEGER IWORK(NGRID,NGRID) C ********** DOMAIN PROCESSOR DECLARATIONS ****** INTEGER N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, A NGRIDX, NGRIDY, LEVEL, QLIMIT,QX(250),QY(250), B GTYPE, PIECE, BGRID, BNEIGH REAL PARAM, EPSGRD, EPSTAN LOGICAL CLOCKW, ARC, FATAL, INHOLE INTEGER TYPE, BPTYPE, HORZ, VERT, BOTH, INTER, JUMP COMMON / DMCINT / N1BND, N1BDPT, NSIDE, IPACKB, MOUTPT, A NGRIDX, NGRIDY, LEVEL, QLIMIT, QX, QY COMMON / DMCREL / PARAM, EPSGRD, EPSTAN COMMON / BNDRYI / IPIECE, NBOUND, NBNDPT COMMON / BNDRYL / CLOCKW, ARC, FATAL, INHOLE COMMON / DMCHAR / TYPE, HORZ, VERT, BOTH, INTER, JUMP C SET PARAMETERS FOR RUN MINPUT = 5 MOUTPT = 6 NGDIMX=NGRID NGDIMY=NGRID NBDIM=300 NPDIM=1 C CREATE GRID DX = 1.0/FLOAT(NGDIMX-1) DY = 1.0/FLOAT(NGDIMY-1) DO 10 I=1,NGDIMX XGRID(I) = FLOAT(I-1)*DX 10 CONTINUE DO 20 I=1,NGDIMY YGRID(I) = FLOAT(I-1)*DY 20 CONTINUE C SET OUTPUT LEVEL LEVEL = 2 C INITIALIZE BOUNDARY FOR OUTER CIRCLE NBOUND = 1 BRANGE(1,1) = 0.0 BRANGE(2,1) = 8.0*ATAN(1.0) CALL REGION(XGRID,YGRID,NGDIMX,NGDIMY,BRANGE,NPDIM,CIRCLE, A .TRUE.,.FALSE.,LEVEL, B GTYPE,XBOUND,YBOUND,PIECE,BPTYPE,BNEIGH, C BGRID,BPARAM,NBDIM,NBPTS,FATAL) C C INITIALIZE BOUNDARY FOR SECOND CIRCLE NBOUND = 2 BRANGE(1,2) = 0.0 BRANGE(2,2) = 8.0*ATAN(1.0) CALL HOLE(XGRID,YGRID,NGDIMX,NGDIMY,BRANGE,NBOUND,CIR2, A .TRUE.,.FALSE.,LEVEL, B GTYPE,XBOUND,YBOUND,PIECE,BPTYPE,BNEIGH, C BGRID,BPARAM,NBDIM,NBPTS,FATAL,IWORK) C C INITIALIZE BOUNDARY FOR THIRD CIRCLE NBOUND = 3 BRANGE(1,3) = 0.0 BRANGE(2,3) = 8.0*ATAN(1.0) CALL HOLE(XGRID,YGRID,NGDIMX,NGDIMY,BRANGE,NBOUND,CIR3, A .TRUE.,.FALSE.,LEVEL, B GTYPE,XBOUND,YBOUND,PIECE,BPTYPE,BNEIGH, C BGRID,BPARAM,NBDIM,NBPTS,FATAL,IWORK) C STOP END C SUBROUTINE CIRCLE(S,X,Y,IPIECE) DATA SMIN, SMAX / 0.0, 6.28319 / IF ( (S.LT.SMIN) .OR. (S.GT.SMAX) ) WRITE(6,10) S 10 FORMAT(5(5H** **),6H S = ,F10.3,14H OUT OF RANGE.) X = 0.5*SIN(S) + 0.5 Y = 0.5*COS(S) + 0.5 RETURN END SUBROUTINE CIR2(S,X,Y,IPIECE) DATA SMIN, SMAX / 0.0, 6.28319 / IF ( (S.LT.SMIN) .OR. (S.GT.SMAX) ) WRITE(6,10) S 10 FORMAT(5(5H** **),6H S = ,F10.3,14H OUT OF RANGE.) X = 0.10*COS(S) + 0.25 Y = 0.10*SIN(S) + 0.25 RETURN END SUBROUTINE CIR3(S,X,Y,IPIECE) DATA SMIN, SMAX / 0.0, 6.28319 / IF ( (S.LT.SMIN) .OR. (S.GT.SMAX) ) WRITE(6,10) S 10 FORMAT(5(5H** **),6H S = ,F10.3,14H OUT OF RANGE.) X = 0.11*COS(S) + 0.68 Y = 0.11*SIN(S) + 0.6 RETURN END