C*********************************************************************** C----PROGRAM CLAF----JULY 18, 1986---- C-----CHANGED ERROR IN SETSYM---- C-----change claf43 to have double precision---- C-----Changed claf41 and claf42 to have different integration C-------for the x and the y axes. C*********************************************************************** PARAMETER (NBIG=180000) DIMENSION P(NBIG) COMMON/CONTRO/NUMFRQ,NLAYER,IFCVAR,IPRNT,IPLOT,IDIM,G,VS,CL COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 COMMON/UNITS/ISCR1,ISCR2,ISCR3 ICLA1=5 ICLA2=6 ICLA3=7 ICLA4=8 ICLA5=9 ISCR1=21 ISCR2=22 OPEN(UNIT=ICLA1,FILE='NEWN1.DAT',STATUS='OLD') OPEN(UNIT=ICLA2,FILE='NEWN2.DAT',STATUS='NEW') OPEN(UNIT=ICLA3,FILE='IL2.DAT',STATUS='OLD',FORM='UNFORMATTED') OPEN(UNIT=ICLA4,FILE='IMPED.DAT',STATUS='NEW',FORM='UNFORMATTED') OPEN(UNIT=ICLA5,FILE='DRFOR.DAT',STATUS='NEW',FORM='UNFORMATTED') C OPEN(UNIT=21,FILE='SCR1.DAT',STATUS='NEW',FORM='UNFORMATTED') OPEN(UNIT=22,FILE='SCR2.DAT',STATUS='NEW',FORM='UNFORMATTED') C*********************************************************************** C-------------------------------------------------READ STATEMENT-------- READ (ICLA1,10) G,VS,CL 10 FORMAT(3E10.3) C-------------------------------------------------READ STATEMENT-------- READ (ICLA1,20) NFDN,NUMFRQ,NCASE,IFCVAR,NGRN,NTYPE,NLAYER 20 FORMAT(16I5) NF6=NFDN*6 WRITE (ICLA2,1001) G,VS,CL WRITE (ICLA2,1003) NFDN,NTYPE,NGRN WRITE (ICLA2,1002) NUMFRQ,NCASE,NLAYER IF(IFCVAR.EQ.1)WRITE (ICLA2,1005) MCNP=1 MCXC=MCNP+NFDN MCYC=MCXC+NFDN MCIS=MCYC+NFDN MCNREP=MCIS+NFDN MCSNX=MCNREP+NFDN MCSNY=MCSNX+NFDN*4 MCKNK=MCSNY+NFDN*4 MCXYBL=MCKNK+NFDN*6 CALL RDFDN1(NFDN,P(MCNP),P(MCXC),P(MCYC),P(MCIS),P(MCNREP), # P(MCSNX),P(MCSNY),P(MCKNK),P(MCXYBL),MEXTRA,SCALE) C----MEXTRA CONTAINS ALL THE ARRAYS OF XB,YB AND LET.---- MCXH=MCXYBL+MEXTRA MCYH=MCXH+NTYPE CALL RDSHP(P(MCXH),P(MCYH),NTYPE,SCALE) MCGF=MCYH+NTYPE C----MAKE MCGF ODD TO ACCOMODATE COMPLEX NUMBERS, AN IBM PROBLEM---- MCGF=(MCGF/2)*2+1 MCK=MCGF+2*4*NGRN CALL MXMEM(NFDN,P(MCNP),P(MCIS),NSTOR,NTR) MCDF=MCK+(NF6**2)*2 MCTR1=MCDF+(NF6*NCASE)*2 MCTR2=MCTR1+(NTR*NF6)*2 ITEMP=NF6 IF(NFDN.EQ.1)ITEMP=0 MCW=MCTR2+(NTR*ITEMP)*2 MWAVES=MCW+NSTOR MCGRD=MWAVES+8*NCASE ITEMP=NCASE IF(NF6.GT.NCASE)ITEMP=NF6 MFINAL=MCGRD+(NTR*ITEMP)*2 ITEMP=MFINAL-1 WRITE (ICLA2,1007) ITEMP 1007 FORMAT(' CORE USAGE=',I10,' SINGLE PRECISION WORDS.') IF(MFINAL.GT.NBIG)GO TO 1010 CALL CLAF2(NFDN,P(MCNP),P(MCXC),P(MCYC), # P(MCIS),P(MCNREP),P(MCSNX), # P(MCSNY),P(MCKNK),P(MCXYBL),P(MCXH),P(MCYH),P(MCGF), # P(MCK),P(MCDF),P(MCTR1),P(MCTR2),P(MCW), # P(MWAVES),P(MCGRD),NTR,NF6,NGRN,NCASE) C=======================================FORMATS FOR OUTPUT============== 1001 FORMAT(4X,28HTHE REFERENCE SHEAR MODULUS=,E11.3// * 4X,34HTHE REFERENCE SHEAR WAVE VELOCITY=,E11.3// * 4X,44HTHE CHARACTERISTIC LENGTH OF THE FOUNDATION=,E11.3/) 1002 FORMAT(4X,7HNUMFRQ=,I5,10X,6HNCASE=,I5,10X,7HNLAYER=,I5/) 1003 FORMAT(' NUMBER OF FOUNDATIONS = ',I2// # ' THE NUMBER OF DIFFERENT SUBREGION SIZES USED = ',I3// # ' THE MAXIMUM SIZE OF THE GREENS FUNCTION TABLE = ',I6/) 1005 FORMAT(4X,46HTHE WAVE VELOCITIES ARE FUNCTIONS OF FREQUENCY/) 1006 FORMAT(4X,33HTHE CORE USAGE OF THIS PROBLEM IS,I7,6H WORDS/ * 4X,29HTHE MAXIMUM CORE ALLOCATED IS,I7,6H WORDS/) C=======================================FORMATS FOR OUTPUT============== 1010 STOP END C*********************************************************************** BLOCK DATA COMMON/PATTRN/ISW(3,4) COMMON/INTEGR/QWI(27),QXI(27),IBN(6) DATA IBN/0,2,5,9,14,20/ DATA QWI/0.50,0.50, * 0.27777778,0.44444444,0.27777778, * 0.17392743,0.32607258,0.32607258,0.17392743, * 0.11846345,0.23931434,0.28444445,0.23931434,0.11846345, * 0.08566225,0.18038079,0.23395697, * 0.23395697,0.18038079,0.08566225, * 0.06474249,0.13985270,0.19091503,0.20897959, * 0.19091503,0.13985270,0.06474249/ DATA QXI/-0.28867514,0.28867514, * -0.38729834, 0.0, 0.38729834, * -0.43056816,-0.16999052, 0.16999052, 0.43056816, * -0.45308993,-0.26923466, 0.0, 0.26923466, 0.45308993, * -0.46623476,-0.33060470,-0.11930960, * 0.11930960, 0.33060470, 0.46623476, * -0.47455396,-0.37076560,-0.20292258, 0.0, * 0.20292258, 0.37076560, 0.47455396/ DATA ISW/1,4,2, 4,1,3, 2,3,1, 3,2,4/ END C*********************************************************************** SUBROUTINE CLAF2(NFDN,NP,XC,YC,IS,NREP,SNX,SNY,KNK,XYBL,XH,YH, # GF,K,DF,TR1,TR2,W,WAVES,GROUND,NTR,NF6,NGRN,NCASE) COMPLEX GF(4,NGRN) COMPLEX K(NF6,NF6),DF(NF6,1),GROUND(NTR,1) DIMENSION NP(NFDN),IS(NFDN),SNX(4,NFDN),SNY(4,NFDN) DIMENSION XC(NFDN),YC(NFDN),KNK(6,NFDN),WAVES(1) DIMENSION XYBL(1),XH(1),YH(1),NREP(NFDN) COMPLEX TR1(NTR,NF6),TR2(NTR,NF6) DIMENSION W(1) COMPLEX G0 COMMON/CONTRO/NUMFRQ,NLAYER,IFCVAR,IPRNT,IPLOT,IDIM,G,VS,CL COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 COMMON/UNITS/ISCR1,ISCR2,ISCR3 COMMON/SGRN/DIV,G0(4),NGRN1 COMMON/SW/ETA PI=3.141592654 c call swbeg ioldt=0 DO 600 NF=1,NUMFRQ DO 20 I=1,NF6 DO 10 J=1,NF6 10 K(I,J)=0.0 DO 20 J=1,NTR 20 TR1(J,I)=0.0 REWIND ISCR2 C-------------------------------------------------READ STATEMENT-------- READ (ICLA1,300) FRQ 300 FORMAT(E10.3) c call swtime(ijklmn) timer1=(ijklmn-ioldt)/60000. ioldt=ijklmn write (*,*) 'Frequency=',frq,timer1 C----------------------------------------------------------------------- C READ IN GREEN'S FUNCTION TABLE. C----------------------------------------------------------------------- IF(NF.NE.1.AND.NLAYER.EQ.1)GO TO 350 IF(NCASE.EQ.0)GO TO 330 IF(NF.NE.1 .AND. IFCVAR.EQ.0)GO TO 330 CALL RDWAVE(WAVES,NCASE) 330 CONTINUE C-------------------------------------------------READ BINARY----------- READ (ICLA3) DIV,NGRN1,FFFR C-------------------------------------------------READ BINARY----------- READ (ICLA3) ((GF(IX,J),IX=1,4),J=1,NGRN1) IF(NGRN1.GT.NGRN)GO TO 700 C----SUBTRACT OUT THE STATIC PART---- DO 340 IX=1,4 G0(IX)=GF(IX,1) DO 340 J=1,NGRN1 340 GF(IX,J)=GF(IX,J)-G0(IX) 350 CONTINUE IF(NLAYER.EQ.1)GO TO 360 RAT=FRQ/FFFR IF(RAT.LT.1.01 .AND. RAT.GT.0.99)GO TO 370 WRITE (ICLA2,1005) FRQ,FFFR GO TO 600 360 IF(FRQ.LE.FFFR)GO TO 370 WRITE (ICLA2,1006) FRQ GO TO 600 370 CONTINUE C----------------------------------------------------------------------- C ACTUAL CALCULATION OF IMPEDANCE MATRIX AND DRIVING FORCES BEGIN. C----------------------------------------------------------------------- ETA=FRQ*2.*PI*CL/VS IF(NFDN.EQ.1)GO TO 800 C----OFF DIAGONAL BLOCKS---- DO 410 IFDN=1,NFDN ISI=IS(IFDN) NPI3=NP(IFDN)*ISI*3 CALL LOCIND(NP,IS,NREP,IFDN,IORD,IXB,IYB,ILET) DO 420 JFDN=IFDN+1,NFDN NPJ3=NP(JFDN)*IS(JFDN)*3 CALL LOCIND(NP,IS,NREP,JFDN,JORD,JXB,JYB,JLET) XCC=XC(IFDN)-XC(JFDN) YCC=YC(IFDN)-YC(JFDN) CALL CLAF4E(W,GF,NPI3,NPJ3,XYBL(IXB),XYBL(IYB), # XYBL(ILET),XYBL(JXB),XYBL(JYB),XYBL(JLET),XH,YH,ISI, # IS(JFDN),SNX(1,IFDN),SNY(1,IFDN),SNX(1,JFDN),SNY(1,JFDN), # XCC,YCC) CALL SAVE(W,NPI3,NPJ3) 420 CONTINUE 410 CONTINUE 800 CONTINUE DO 810 IFDN=1,NFDN ISI=IS(IFDN) NPPI3=NP(IFDN)*3 NPI3=NP(IFDN)*ISI*3 MEG=1 METMP=MEG+NPI3*NPPI3*2 MEORD=METMP+NPI3*2 MERHS=MEORD+NPI3 CALL LOCIND(NP,IS,NREP,IFDN,IORD,IXB,IYB,ILET) IF(NREP(IFDN).EQ.1)GO TO 820 CALL CLAF4D(W(MEG),W(MERHS),GF,NPPI3,XYBL(IXB),XYBL(IYB), # XYBL(ILET),XH,YH,ISI,SNX(1,IFDN),SNY(1,IFDN)) CALL DECOMP(W(MEG),W(METMP),W(MEORD),NPPI3,ISI,SNX(1,IFDN)) CALL SAVE(W(MEG),NPPI3,NPI3) CALL SAVEI(W(MEORD),NPI3) CALL ZEROTH(W(MEG),W(MERHS),W(MEORD),W(METMP),NPPI3, # ISI,KNK(1,IFDN)) 820 CALL KEEPTR(W(MERHS),TR1,IORD,NPPI3,NTR,NF6,IFDN,ISI, # SNX(1,IFDN),SNY(1,IFDN),KNK(1,IFDN)) 810 CONTINUE IF(NFDN.EQ.1)GO TO 531 C----ITERATION BEGINS---- CALL COPYAR(TR1,GROUND,(NTR*NF6)) DIFF=1.0 DO 500 ITER=1,20 REWIND ISCR2 CALL MULT1(TR1,TR2,W,NTR,NF6,NFDN,NP,IS,NREP,ITER) DO 550 IFDN=1,NFDN ISI=IS(IFDN) NPPI3=NP(IFDN)*3 NPI3=NPPI3*ISI CALL LOCIND(NP,IS,NREP,IFDN,IORD,IXB,IYB,ILET) MEG=1 METMP=MEG+NPI3*NPPI3*2 MEORD=METMP+NPI3*2 MERHS=MEORD+NPI3 CALL SOLVE(W(MEG),TR2,W(MEORD),W(METMP),IORD,NREP(IFDN), # ISI,SNX(1,IFDN),SNY(1,IFDN),NPI3,NF6,NTR) 550 CONTINUE CALL ADDAR(TR2,GROUND,(NTR*NF6)) IF(ITER.LT.3)GO TO 546 CALL COMPAR(TR2,GROUND,(NTR*NF6),DIFF) C WRITE (ICLA2,1050) ITER,DIFF 1050 FORMAT(' ITERATION NUMBER ',I2,' AVERAGE ERROR=',F8.5) IF(DIFF.LT.0.02)GO TO 530 546 CALL COPYAR(TR2,TR1,(NTR*NF6)) 500 CONTINUE 530 CALL COPYAR(GROUND,TR1,(NTR*NF6)) 531 DO 540 IFDN=1,NFDN ISI=IS(IFDN) NPPI3=NP(IFDN)*3 NPI3=NPPI3*ISI CALL LOCIND(NP,IS,NREP,IFDN,IORD,IXB,IYB,ILET) CALL FORMK(K,TR1,IFDN,2,IORD,NTR,NF6,NPI3,XYBL(IXB),XYBL(IYB), # XYBL(ILET),XH,YH,ISI,SNX(1,IFDN),SNY(1,IFDN)) IF(NCASE.EQ.0)GO TO 540 CALL CLAF3(GROUND,WAVES,NTR,NCASE,IORD,NPI3, # XYBL(IXB),XYBL(IYB),XYBL(ILET),XH,YH,ISI, # SNX(1,IFDN),SNY(1,IFDN),XC(IFDN),YC(IFDN)) 540 CONTINUE CALL WRK(K,NF6,ETA) IF(NCASE.EQ.0)GO TO 600 CALL FORMDF(TR1,GROUND,DF,NTR,NF6,NCASE) CALL WRDF(DF,NF6,NCASE,ETA) 600 CONTINUE C=======================================FORMATS FOR OUTPUT============== 1003 FORMAT(/2X,25HINCIDENT WAVE INPUT DATA:,// * 4X,6HBETA/C,5X,5HTHETA,3X,12HLONGITUDINAL,5X, * 10HTRANSVERSE,7X,8HVERTICAL) 1004 FORMAT(F10.4,F10.2,6F8.3) 1005 FORMAT(4X,13HTHE FREQUENCY,E11.3,24H DOES NOT MATCH THAT OF , * 9HTHE TABLE,E11.3,42H THE FREQUENCY INPUT MAY BE OUT OF ORDER.) 1006 FORMAT(4X,13HTHE FREQUENCY,E11.3,29H EXCEEDS THAT OF THE MAXIMUM, * 51HFREQUENCY, GREENS FUNCTION TABLE IS NOT SUFFICIENT.) 1007 FORMAT(F7.2,I5,5X,48HIMPEDANCE MATRIX FOLLOWS IN 3(2X,2E10.3) FORM *AT.) 1008 FORMAT(3(2X,2E10.3)) 1009 FORMAT(F7.2,I5,5X,43HINPUT MOTIONS FOLLOW IN 3(2X,2E10.3) FORMAT) 1010 FORMAT(1H1,5X,26HSUMMARY OF IMPORTANT TERMS,/) 1011 FORMAT(/5X,27HSTIFFNESS TERMS FOR D.O.F.S,6I3) 1012 FORMAT(2X,3HA0=,F7.3,3X,6E10.3) 1013 FORMAT(/5X,25HDAMPING TERMS FOR D.O.F.S,6I3) 1014 FORMAT(/5X,21HINCIDENT WAVE NUMBER:,I5/ * 5X,38HAMPLITUDES OF INPUT MOTION FOR D.O.F.S,6I3) C=======================================FORMATS FOR OUTPUT============== 700 RETURN END C*********************************************************************** SUBROUTINE CLAF3(RHS,WAVES,NTR,NCASE,IORD, # NP3,XB,YB,LET,XH,YH,IS,SNX,SNY,XC,YC) COMPLEX RHS(NTR,NCASE),U(3) COMMON/SW/ETA DIMENSION SNX(4),SNY(4),WAVES(1) DIMENSION XB(1),YB(1),XH(1),YH(1),LET(1) NP=NP3/(3*IS) DO 300 KSI=1,IS DO 300 IM=1,NP I=(IM-1)*3 + (KSI-1)*NP*3 +(IORD-1) XL=XB(IM)*SNX(KSI) + XC YL=YB(IM)*SNY(KSI) + YC IG=LET(IM) HX=XH(IG) HY=YH(IG) DO 100 IW=1,NCASE INDEX=(IW-1)*8+1 CALL CLAF31(XL,YL,HX,HY,WAVES(INDEX),U) DO 50 IR=1,3 50 RHS(IR+I,IW)=U(IR) 100 CONTINUE 300 CONTINUE RETURN END C*********************************************************************** SUBROUTINE CLAF31(X,Y,HX,HY,WAVES,U) COMMON/SW/ETA COMMON/INTEGR/QWI(27),QXI(27),IBN(6) COMPLEX U(3),UL,UT,UV DIMENSION WAVES(8) C======================================================================= C CLAF31 CALCULATES THE AVERAGE OF THE FREE FIELD MOTION IN A C SUBREGION. THE WAVE FORMS ARE CURRENTLY ASSUMED TO BE PLANE WAVES, C IT CAN, HOWEVER, BE REPLACED BY OTHER WAVE FORMS IF THE USER C SUPPLIES A NEW SUBROUTINE. C======================================================================= BOC=WAVES(1) TH=WAVES(2) UL=CMPLX(WAVES(3),WAVES(4)) UT=CMPLX(WAVES(5),WAVES(6)) UV=CMPLX(WAVES(7),WAVES(8)) C----TAKE AVERAGE OF PLANE WAVE FUNCTION---- U(3)=0.0 DO 10 I=1,2 XX=X+QXI(I)*HX DO 10 J=1,2 YY=Y+QXI(J)*HY RI=ETA*BOC*(XX*COS(TH)+YY*SIN(TH)) U(3)=U(3)+COS(RI)-(0.,1.)*SIN(RI) 10 CONTINUE U(3)=U(3)*HX*HY/4. U(1)=(UL*COS(TH)-UT*SIN(TH))*U(3) U(2)=(UL*SIN(TH)+UT*COS(TH))*U(3) U(3)=UV*U(3) RETURN END C*********************************************************************** SUBROUTINE CLAF4D(A,ALPHA,GF,NPP3,XB,YB,LET,XH,YH,IS,SNX,SNY) COMMON/SW/ETA COMPLEX GF(4,1) COMPLEX UH(6),A(NPP3,1),ALPHA(NPP3,6) DIMENSION SNX(4),SNY(4) DIMENSION XB(1),YB(1),XH(1),YH(1),LET(1) C======================================================================= C CLAF4 SETS UP THE INFLUENCE MATRIX, IT CALLS CLAF41 TO CALCULATE C THE INFLUENCE FUNCTIONS. C======================================================================= NP=NPP3/3 DO 300 IM=1,NP I=(IM-1)*3 I1=I+1 I2=I+2 I3=I+3 XL=XB(IM) YL=YB(IM) IG=LET(IM) HXO=XH(IG) HYO=YH(IG) C----SET UP RIGHT HAND SIDE VECTOR---- ZAREA=HXO*HYO DO 110 JRHS=1,6 ALPHA(I1,JRHS)=0.0 ALPHA(I2,JRHS)=0.0 110 ALPHA(I3,JRHS)=0.0 ALPHA(I1,1)=ZAREA ALPHA(I2,2)=ZAREA ALPHA(I3,3)=ZAREA ALPHA(I1,6)=-YL*ZAREA ALPHA(I2,6)= XL*ZAREA ALPHA(I3,4)= YL*ZAREA ALPHA(I3,5)=-XL*ZAREA DO 150 KSJ=1,IS IST=1 IF(KSJ.EQ.1)IST=IM DO 150 JM=IST,NP J=(JM-1)*3 + (KSJ-1)*NP*3 J1=J+1 J2=J+2 J3=J+3 XO=XL-XB(JM)*SNX(KSJ) YO=YL-YB(JM)*SNY(KSJ) IG=LET(JM) HX=XH(IG) HY=YH(IG) CALL CLAF41(ETA,XO,YO,HX,HY,HXO,HYO,UH,GF) A(I1,J1)= UH(1) A(I1,J2)= UH(2) A(I1,J3)= UH(3) A(I2,J1)= UH(2) A(I2,J2)= UH(4) A(I2,J3)= UH(5) A(I3,J1)=-UH(3) A(I3,J2)=-UH(5) A(I3,J3)= UH(6) 150 CONTINUE 300 CONTINUE C--DO REFLECTION---- DO 400 I=1,NPP3-1 DO 400 J=I+1,NPP3 400 A(J,I)=A(I,J) RETURN END C*********************************************************************** SUBROUTINE CLAF4E(A,GF,NP13,NP23,XB1,YB1,LET1,XB2,YB2,LET2,XH,YH, # IS1,IS2,SNX1,SNY1,SNX2,SNY2,XCC,YCC) COMMON/SW/ETA COMPLEX GF(4,1) COMPLEX UH(6),A(NP13,NP23) DIMENSION XB1(1),YB1(1),XB2(1),YB2(1),XH(1),YH(1),LET1(1),LET2(1) DIMENSION SNX1(4),SNY1(4),SNX2(4),SNY2(4) C======================================================================= C CLAF4 SETS UP THE INFLUENCE MATRIX, IT CALLS CLAF41 TO CALCULATE C THE INFLUENCE FUNCTIONS. C======================================================================= NP1=NP13/(3*IS1) NP2=NP23/(3*IS2) DO 300 KSI=1,IS1 DO 300 IM=1,NP1 I=(IM-1)*3 + (KSI-1)*NP1*3 XL=XB1(IM)*SNX1(KSI) + XCC YL=YB1(IM)*SNY1(KSI) + YCC IG=LET1(IM) HXO=XH(IG) HYO=YH(IG) DO 150 KSJ=1,IS2 DO 150 JM=1,NP2 J=(JM-1)*3 + (KSJ-1)*NP2*3 XO=XL-XB2(JM)*SNX2(KSJ) YO=YL-YB2(JM)*SNY2(KSJ) IG=LET2(JM) HX=XH(IG) HY=YH(IG) CALL CLAF41(ETA,XO,YO,HX,HY,HXO,HYO,UH,GF) I1=I+1 I2=I+2 I3=I+3 J1=J+1 J2=J+2 J3=J+3 A(I1,J1)= UH(1) A(I1,J2)= UH(2) A(I1,J3)= UH(3) A(I2,J1)= UH(2) A(I2,J2)= UH(4) A(I2,J3)= UH(5) A(I3,J1)=-UH(3) A(I3,J2)=-UH(5) A(I3,J3)= UH(6) 150 CONTINUE 300 CONTINUE RETURN END C*********************************************************************** SUBROUTINE CLAF41(ETA,XO,YO,HX,HY,HXO,HYO,UH,GF) COMPLEX GF(4,1),UH(6),WH(6) COMMON/INTEGR/QWI(27),QXI(27),IBN(6) C======================================================================= C CLAF41 CALCULATES THE INFLUENCE FUNCTIONS. A 4-POINT GAUSSIAN C QUADRATURE IS USED TO CALCULATE THE AVERAGE. C======================================================================= XTEST=ETA*HX YTEST=ETA*HY NUMX=XTEST NUMX=MAX0(2,NUMX) NUMY=YTEST NUMY=MAX0(2,NUMY) IF(XO.LT.HX .AND. YO.LT.HY)THEN NUMX=6 NUMY=6 ENDIF NUMX=MIN0(7,NUMX) NUMY=MIN0(7,NUMY) DO 110 JJ=1,6 110 UH(JJ)=0. DO 130 I=1,2 XZ=XO+QXI(I)*HXO DO 130 J=1,2 YZ=YO+QXI(J)*HYO CALL CLAF42(ETA,XZ,YZ,HX,HY,WH,GF,NUMX,NUMY) DO 120 JJ=1,6 120 UH(JJ)=UH(JJ)+WH(JJ) 130 CONTINUE ETEST=HXO*HYO/4. DO 140 JJ=1,6 140 UH(JJ)=UH(JJ)*ETEST RETURN END C*********************************************************************** SUBROUTINE CLAF42(ETA,X0,Y0,HHX,HHY,P,GF,NUMX,NUMY) COMPLEX P(6),U1,U2,U3,U4,GF(4,1),G0 COMMON/SGRN/DIV,G0(4),NGRN COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 COMMON/INTEGR/QWI(27),QXI(27),IBN(6) COMMON/FASTAC/XQ,YQ,XQ2,YQ2,D,D2,D3,WW C======================================================================= C CLAF42 CALCULATES THE THREE DISPLACEMENTS AT A POINT DUE TO C UNIFORM LOADS AT A RECTANGULAR SUBREGION. C======================================================================= C P(1)=G11 P(2)=G12 P(3)=G13 C P(4)=G22 P(5)=G23 C P(6)=G33 CALL CLAF43(X0,Y0,HHX/2.,HHY/2.,P) IF(ETA.EQ.0.0)RETURN JBX=IBN(NUMX-1) JBY=IBN(NUMY-1) DO 30 IY=1,NUMY JBCY=JBY+IY YQ=Y0+QXI(JBCY)*HHY WW1=QWI(JBCY)*HHX*HHY YQ2=YQ*YQ DO 20 IX=1,NUMX JBCX=JBX+IX XQ=X0+QXI(JBCX)*HHX WW=WW1*QWI(JBCX) XQ2=XQ*XQ D2=XQ2+YQ2 IF(D2.LT.0.00001)GO TO 10 D=SQRT(D2) D3=D*D2 A=ETA*D N1=A/DIV+1 IF((N1+1).GT.NGRN)GO TO 40 AP=(A-(N1-1)*DIV)/DIV CALL CLAF44(AP,GF(1,N1),P) 10 CONTINUE 20 CONTINUE 30 CONTINUE RETURN 40 WRITE (ICLA2,1000) 1000 FORMAT(2X,41HGREEN'S FUNCTION TABLE HAS BEEN EXCEEDED.) STOP END C*********************************************************************** SUBROUTINE CLAF43(X0,Y0,B,C,US) COMPLEX US(6),G0 DOUBLE PRECISION A1,A2,A3,A4,T13,T23,T14,T24 COMMON/SGRN/DIV,G0(4),NGRN C======================================================================= C CLAF43 CALCULATES THE STATIC PART OF THE INFLUENCE FUNCTIONS. C======================================================================= A1=B-X0 A2=-B-X0 A3=C-Y0 A4=-C-Y0 IA1=A1*10000. IA2=A2*10000. IA3=A3*10000. IA4=A4*10000. T13=DSQRT(A1**2+A3**2) T23=DSQRT(A2**2+A3**2) T14=DSQRT(A1**2+A4**2) T24=DSQRT(A2**2+A4**2) Q32=0. Q42=0. Q13=0. Q23=0. Q15=0. Q25=0. Q35=0. Q45=0. Q16=0. Q26=0. Q36=0. Q46=0. 100 IF(IA1.EQ.0)GO TO 200 Q13=A1*DLOG((A3+T13)/(A4+T14)) Q15=2.*A1*(DATAN(A3/A1)-DATAN(A4/A1)) Q16=2.*A1*DLOG(T13/T14) 200 IF(IA2.EQ.0)GO TO 300 Q23=-A2*DLOG((A3+T23)/(A4+T24)) Q25=-2.*A2*(DATAN(A3/A2)-DATAN(A4/A2)) Q26=-2.*A2*DLOG(T23/T24) 300 IF(IA3.EQ.0)GO TO 400 Q32=A3*DLOG((A1+T13)/(A2+T23)) Q35=2.*A3*DLOG(T13/T23) Q36=2.*A3*(DATAN(A1/A3)-DATAN(A2/A3)) 400 IF(IA4.EQ.0)GO TO 500 Q42=-A4*DLOG((A1+T14)/(A2+T24)) Q45=-2.*A4*DLOG(T14/T24) Q46=-2.*A4*(DATAN(A1/A4)-DATAN(A2/A4)) 500 CONTINUE S2=Q32+Q42 S3=Q13+Q23 S1=S2+S3 S4=-T13+T23+T14-T24 S5=-0.5*(Q15+Q25+Q35+Q45) S6=-0.5*(Q16+Q26+Q36+Q46) US(1)=G0(1)*S2-G0(2)*S3 US(2)=(G0(1)+G0(2))*S4 US(3)=-G0(3)*S5 US(4)=G0(1)*S3-G0(2)*S2 US(5)=-G0(3)*S6 US(6)=G0(4)*S1 RETURN END C*********************************************************************** SUBROUTINE CLAF44(AP,GF,P) COMMON/FASTAC/XQ,YQ,XQ2,YQ2,D,D2,D3,WW COMPLEX P(6),GF(8),U1,U2,U3,U4 F1=(1.-AP) F2=AP U1=F1*GF(1)+F2*GF(5) U2=F1*GF(2)+F2*GF(6) U3=F1*GF(3)+F2*GF(7) U4=F1*GF(4)+F2*GF(8) P(1)=P(1)+WW*(U1*XQ2-U2*YQ2)/D3 P(2)=P(2)+WW*(U1+U2)*XQ*YQ/D3 P(3)=P(3)-WW*U3*XQ/D2 P(4)=P(4)+WW*(U1*YQ2-U2*XQ2)/D3 P(5)=P(5)-WW*U3*YQ/D2 P(6)=P(6)+WW*U4/D RETURN END C*********************************************************************** SUBROUTINE RDFDN1(NFDN,NP,XC,YC,IS,NREP,SNX,SNY,KNK,A, # MEXTRA,SCALE) COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 LOGICAL CHECK DIMENSION NP(NFDN),IS(NFDN),SNX(4,NFDN),SNY(4,NFDN),A(1) DIMENSION XC(NFDN),YC(NFDN),KNK(6,NFDN),NREP(NFDN) MEXTRA=0 READ (ICLA1,1) SCALE 1 FORMAT(F10.5) WRITE (ICLA2,1002) SCALE DO 10 IFDN=1,NFDN WRITE (ICLA2,1001) IFDN READ (ICLA1,2) NPI,IFSX,IFSY,NREP(IFDN) 2 FORMAT(4I5) READ (ICLA1,3) XC(IFDN),YC(IFDN) 3 FORMAT(2F10.5) XC(IFDN)=XC(IFDN)*SCALE YC(IFDN)=YC(IFDN)*SCALE WRITE (ICLA2,1003) XC(IFDN),YC(IFDN) NP(IFDN)=NPI CALL SETSYM(IFSX,IFSY,IS(IFDN),SNX(1,IFDN), # SNY(1,IFDN),KNK(1,IFDN)) NPI3=NP(IFDN)*IS(IFDN) WRITE (ICLA2,1004) NPI3 IF(NREP(IFDN).EQ.1)GO TO 20 MDXB=MEXTRA+1 MDYB=MDXB+NPI MDLET=MDYB+NPI MEXTRA=MEXTRA+3*NPI CALL RDFDN2(A(MDXB),A(MDYB),A(MDLET),NPI,SCALE) GO TO 10 20 IF(IFDN.NE.1)GO TO 30 WRITE (ICLA2,1020) STOP 30 WRITE (ICLA2,1005) IF(NPI.NE.NP(IFDN-1))WRITE (ICLA2,1021) IFDN 10 CONTINUE 1002 FORMAT(/' FACTOR USED TO SCALE INPUT DATA = ',F10.5) 1001 FORMAT(/' >>>FOUNDATION NUMBER ',I2) 1003 FORMAT(' LOCATION OF THE ORIGIN = (',F10.5,',',F10.5,')'/) 1004 FORMAT(/' NUMBER OF SUBREGIONS = ',I4) 1005 FORMAT(' THE GEOMETRY OF THIS FOUNDATION IS THE SAME AS ', # 'THE PREVIOUS FOUNDATION.') 1020 FORMAT(' CANNOT HAVE NREP=1 FOR THE FIRST FOUNDATION.') 1021 FORMAT( ' THE NREP VALUE OF FOUNDATION NUMBER',I5, # ' INDICATES THAT IT SHOULD BE A REPLICA OF THE PREVIOUS ', # 'FOUNDATION, BUT THE VALUES OF NP DO NOT MATCH.') RETURN END C*********************************************************************** SUBROUTINE RDFDN2(XB,YB,LET,NP,SCALE) COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 DIMENSION XB(NP),YB(NP),LET(NP) DO 10 I=1,NP READ (ICLA1,1) XB(I),YB(I),LET(I) 1 FORMAT(2F10.5,I5) XB(I)=XB(I)/SCALE YB(I)=YB(I)/SCALE WRITE (ICLA2,1001) I,XB(I),YB(I),LET(I) 10 CONTINUE 1001 FORMAT(3X,11HSUBREGION #,I3,10H:CENTROID,,2F7.3,6H, TYPE,I3) RETURN END C*********************************************************************** SUBROUTINE RDSHP(XH,YH,NTYPE,SCALE) COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 DIMENSION XH(NTYPE),YH(NTYPE) WRITE (ICLA2,1001) DO 10 I=1,NTYPE READ (ICLA1,1) XH(I),YH(I) 1 FORMAT(2F10.5) XH(I)=XH(I)/SCALE YH(I)=YH(I)/SCALE WRITE (ICLA2,1002) I,XH(I),YH(I) 10 CONTINUE 1001 FORMAT(/' DEFINITION OF SUBREGION DIMENSIONS:'/) 1002 FORMAT(3X,28HDIMENSIONS OF SUBREGION TYPE,I3,3H IS,2F7.3) RETURN END C*********************************************************************** SUBROUTINE MXMEM(NFDN,NP,IS,NSTOR,NTR) DIMENSION NP(NFDN),IS(NFDN) NLMAX=0 NTR=0 DO 10 IFDN=1,NFDN NPT1=NP(IFDN)*3 NPT2=NPT1*IS(IFDN) NPS=2*NPT1*NPT2 + 2*NPT2 + NPT2 + 2*(NPT1*6) IF(NPS.GT.NLMAX)NLMAX=NPS NTR=NTR+NPT2 DO 5 JFDN=1,NFDN IF(IFDN.EQ.JFDN)GO TO 5 NPQ=NPT2*(NP(JFDN)*IS(JFDN)*3)*2 IF(NPQ.GT.NLMAX)NLMAX=NPQ 5 CONTINUE 10 CONTINUE NSTOR=NLMAX RETURN END C*********************************************************************** SUBROUTINE LOCIND(NP,IS,NREP,IFDN,IORD,IXB,IYB,ILET) DIMENSION NP(1),IS(1),NREP(1) IORD=1 IF(IFDN.EQ.1)GO TO 20 DO 10 I=1,IFDN-1 10 IORD=IORD+NP(I)*IS(I)*3 20 ITMP=1 IF(IFDN.EQ.1)GO TO 40 DO 30 I=1,IFDN-1 IF(NREP(I).EQ.1)GO TO 30 ITMP=ITMP+NP(I)*3 30 CONTINUE IF(NREP(IFDN).NE.1)GO TO 40 ITMP=ITMP-NP(IFDN-1)*3 40 IXB=ITMP IYB=IXB+NP(IFDN) ILET=IYB+NP(IFDN) RETURN END C*********************************************************************** SUBROUTINE ZEROTH(G,RHS,ORDER,TMP,NPP3,IS,KNK) DIMENSION G(1),TMP(1),KNK(6) COMPLEX RHS(NPP3,6) INTEGER ORDER(1) IF(IS.GT.1)GO TO 100 CALL LUSOLV(G,NPP3,RHS,NPP3,6,1,ORDER,TMP) RETURN 100 DO 200 IRIG=1,6 ISET=KNK(IRIG) MEG=(ISET-1)*(2*NPP3**2)+1 MEORD=(ISET-1)*NPP3+1 CALL LUSOLV(G(MEG),NPP3,RHS(1,IRIG),NPP3,1,1, # ORDER(MEORD),TMP) 200 CONTINUE RETURN END C*********************************************************************** SUBROUTINE KEEPTR(RHS,TR1,IORD,NPP3,NTR,NF6,IFDN, # IS,SNX,SNY,KNK) COMMON/PATTRN/ISW(3,4) COMPLEX RHS(NPP3,6),TR1(NTR,NF6) DIMENSION S(4),SNX(4),SNY(4),KNK(6) II=IORD-1 JJ=(IFDN-1)*6 S(1)=1.0 DO 40 J=1,6 ISET=KNK(J) DO 30 ICOMP=1,3 IPAT=ISW(ICOMP,ISET) S(2)=SNX(IPAT) S(3)=SNY(IPAT) S(4)=S(2)*S(3) DO 20 IP=ICOMP,NPP3,3 DO 10 KQUAD=1,IS KQ1=(KQUAD-1)*NPP3+II 10 TR1(KQ1+IP,JJ+J)=RHS(IP,J)*S(KQUAD) 20 CONTINUE 30 CONTINUE 40 CONTINUE RETURN END C*********************************************************************** SUBROUTINE FORMK(K,TR,IFDN,ITER,IORD,NTR,NF6,NP3,XB,YB,LET,XH, # YH,IS,SNX,SNY) COMPLEX K(NF6,NF6),TR(NTR,NF6) DIMENSION XB(1),YB(1),LET(1),XH(1),YH(1),SNX(4),SNY(4) NP=NP3/(3*IS) IF1=(IFDN-1)*6+1 IF2=IF1+1 IF3=IF2+1 IF4=IF3+1 IF5=IF4+1 IF6=IF5+1 JB=1 JE=NF6 IF(ITER.EQ.1)THEN JB=IF1 JE=IF6 ENDIF DO 100 KSI=1,IS DO 50 IM=1,NP I1=(IM-1)*3 + (KSI-1)*NP*3 + IORD I2=I1+1 I3=I2+1 XL=XB(IM)*SNX(KSI) YL=YB(IM)*SNY(KSI) IG=LET(IM) ZAREA=XH(IG)*YH(IG) DO 30 J=JB,JE K(IF1,J)=K(IF1,J)+ZAREA*TR(I1,J) K(IF2,J)=K(IF2,J)+ZAREA*TR(I2,J) K(IF3,J)=K(IF3,J)+ZAREA*TR(I3,J) K(IF4,J)=K(IF4,J)+YL*ZAREA*TR(I3,J) K(IF5,J)=K(IF5,J)-XL*ZAREA*TR(I3,J) K(IF6,J)=K(IF6,J)+(-YL*TR(I1,J)+XL*TR(I2,J))*ZAREA 30 CONTINUE 50 CONTINUE 100 CONTINUE RETURN END C*********************************************************************** SUBROUTINE MULT1(TR1,TR2,G,NTR,NF6,NFDN,NP,IS,NREP,ITER) COMMON/UNITS/ISCR1,ISCR2,ISCR3 COMPLEX TR1(NTR,NF6),TR2(NTR,NF6),G(1) DIMENSION NP(NFDN),IS(NFDN),NREP(NFDN) DO 100 I=1,NTR DO 100 J=1,NF6 100 TR2(I,J)=0.0 DO 1000 IFDN=1,NFDN-1 CALL LOCIND(NP,IS,NREP,IFDN,IORD,ITM1,ITM2,ITM3) NROW=NP(IFDN)*IS(IFDN)*3 DO 500 JFDN=IFDN+1,NFDN CALL LOCIND(NP,IS,NREP,JFDN,JORD,ITM1,ITM2,ITM3) NCOL=NP(JFDN)*IS(JFDN)*3 CALL RECALU(G,NCOL,NROW) CALL MULT2(TR1,TR2,G,NTR,NF6,IORD,JORD,NCOL,NROW, # ITER,JFDN) 500 CONTINUE 1000 CONTINUE REWIND ISCR2 DO 2000 JFDN=1,NFDN-1 CALL LOCIND(NP,IS,NREP,JFDN,JORD,ITM1,ITM2,ITM3) NCOL=NP(JFDN)*IS(JFDN)*3 DO 1500 IFDN=JFDN+1,NFDN CALL LOCIND(NP,IS,NREP,IFDN,IORD,ITM1,ITM2,ITM3) NROW=NP(IFDN)*IS(IFDN)*3 CALL RECALL(G,NCOL,NROW) CALL MULT2(TR1,TR2,G,NTR,NF6,IORD,JORD,NCOL,NROW, # ITER,JFDN) 1500 CONTINUE 2000 CONTINUE RETURN END C*********************************************************************** SUBROUTINE MULT2(TR1,TR2,G,NTR,NF6,IORD,JORD,NCOL,NROW, # ITER,JFDN) COMPLEX TR1(NTR,NF6),TR2(NTR,NF6),G(NCOL,NROW),SUM II=(IORD-1) KB=1 KE=NF6 IF(ITER.EQ.1)THEN KE=JFDN*6 KB=KE-6+1 ENDIF DO 100 I=1,NROW DO 100 K=KB,KE SUM=TR2(II+I,K) CALL MULT3(G(1,I),TR1(JORD,K),NCOL,SUM) TR2(II+I,K)=SUM 100 CONTINUE RETURN END C*********************************************************************** SUBROUTINE MULT3(G,TR1,NN,SUM) COMPLEX G(NN),TR1(NN),SUM DO 1 J=1,NN 1 SUM=SUM-G(J)*TR1(J) RETURN END C*********************************************************************** SUBROUTINE WRK(K,NF6,ETA) COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 COMPLEX K(NF6,NF6) WRITE (ICLA4) ETA,NF6,NF6 DO 10 I=1,NF6 WRITE (ICLA4) (K(I,J),J=1,NF6) 10 CONTINUE RETURN END C*********************************************************************** SUBROUTINE WRDF(DF,NF6,NCASE,ETA) COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 COMPLEX DF(NF6,NCASE) WRITE (ICLA5) ETA,NF6,NCASE DO 10 I=1,NF6 WRITE (ICLA5) (DF(I,J),J=1,NCASE) 10 CONTINUE RETURN END C*********************************************************************** SUBROUTINE COPYAR(A,B,N) COMPLEX A(N),B(N) DO 1 I=1,N 1 B(I)=A(I) RETURN END C*********************************************************************** SUBROUTINE ADDAR(A,B,N) COMPLEX A(N),B(N) DO 1 I=1,N 1 B(I)=B(I)+A(I) RETURN END C*********************************************************************** SUBROUTINE COMPAR(TR2,TR3,N,SUM) COMPLEX TR2(N),TR3(N) SUM=0.0 DO 1 I=1,N 1 SUM=SUM+CABS(TR2(I)/TR3(I)) SUM=SUM/FLOAT(N) RETURN END C*********************************************************************** SUBROUTINE FORMDF(TR1,RHS,DF,NTR,NF6,NCASE) COMPLEX DF(NF6,NCASE),TR1(NTR,NF6),RHS(NTR,NCASE) COMPLEX SUM DO 10 I=1,NF6 DO 10 J=1,NCASE SUM=0.0 DO 5 KK=1,NTR 5 SUM=SUM+TR1(KK,I)*RHS(KK,J) DF(I,J)=SUM 10 CONTINUE RETURN END C*********************************************************************** SUBROUTINE RDWAVE(W,NCASE) COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 DIMENSION W(1) INDEX=0 DO 10 I=1,NCASE READ (ICLA1,1) BOC,TH,ULR,ULI,UTR,UTI,UVR,UVI 1 FORMAT(F10.4,F10.2,6F8.3) W(INDEX+1)=BOC W(INDEX+2)=TH*3.141592654/180. W(INDEX+3)=ULR W(INDEX+4)=ULI W(INDEX+5)=UTR W(INDEX+6)=UTI W(INDEX+7)=UVR W(INDEX+8)=UVI INDEX=INDEX+8 10 CONTINUE RETURN END C*********************************************************************** SUBROUTINE DECOMP(G,TMP,ORDER,NPPI3,IS,SNX) DIMENSION G(1),TMP(1),SNX(4) INTEGER ORDER(1) IF(IS.GT.1) # CALL FOLDL(G,NPPI3,IS,SNX) MEORD1=1 MEG1=1 DO 10 KSI=1,IS CALL LUDCMP(G(MEG1),NPPI3,ORDER(MEORD1),TMP) MEG1=MEG1+(NPPI3**2)*2 MEORD1=MEORD1+NPPI3 10 CONTINUE RETURN END C*********************************************************************** SUBROUTINE SOLVE(G,TR2,ORDER,TMP,IORD,NREPI, # IS,SNX,SNY,NPI3,NF6,NTR) DIMENSION G(1),TMP(1),TR2(NTR,NF6) DIMENSION SNX(4),SNY(4) INTEGER ORDER(1) NPPI3=NPI3/IS IF(NREPI.NE.1)THEN CALL RECALL(G,NPPI3,NPI3) CALL RECALI(ORDER,NPI3) ENDIF IF(IS.GT.1) # CALL FOLDR(TR2,NPI3,NTR,NF6,IS,SNX,IORD) MEORD1=1 MEG1=1 KROWS=IORD DO 10 KSI=1,IS CALL LUSOLV(G(MEG1),NPPI3,TR2,NTR,NF6,KROWS, # ORDER(MEORD1),TMP) KROWS=KROWS+NPPI3 MEORD1=MEORD1+NPPI3 MEG1=MEG1+(NPPI3**2)*2 10 CONTINUE IF(IS.GT.1) # CALL UNFOLD(TR2,NPI3,NTR,NF6,TMP,IS,SNX,SNY,IORD) RETURN END C*********************************************************************** SUBROUTINE SAVE(R,NROW,NCOL) COMMON/UNITS/ISCR1,ISCR2,ISCR3 COMPLEX R(NROW,NCOL) C----WRITE OUT BY COLUMN---- WRITE (ISCR2) R RETURN END C*********************************************************************** SUBROUTINE SAVEI(INTE,NROW) COMMON/UNITS/ISCR1,ISCR2,ISCR3 INTEGER INTE(NROW) WRITE (ISCR2) INTE RETURN END C*********************************************************************** SUBROUTINE RECALL(GR,NROW,NCOL) COMMON/UNITS/ISCR1,ISCR2,ISCR3 C----READ IN BY ROW, THE SAME WAY IT WAS WRITTEN---- COMPLEX GR(NROW,NCOL) READ (ISCR2) GR RETURN END C*********************************************************************** SUBROUTINE RECALU(GR,NCOL,NROW) COMMON/UNITS/ISCR1,ISCR2,ISCR3 C----READ IN AS THE TRANSPOSE OF THE WRITTEN MATRIX---- COMPLEX GR(NCOL,NROW) READ (ISCR2) ((GR(J,I),I=1,NROW),J=1,NCOL) RETURN END C*********************************************************************** SUBROUTINE RECALI(INTE,NROW) COMMON/UNITS/ISCR1,ISCR2,ISCR3 INTEGER INTE(NROW) READ (ISCR2) INTE RETURN END C*********************************************************************** SUBROUTINE SETSYM(IFSX,IFSY,IS,SNX,SNY,KNK) COMMON/UNITC/ICLA1,ICLA2,ICLA3,ICLA4,ICLA5 DIMENSION SNX(4),SNY(4),KNK(6) IS=1 C----------------------------------------------------------------------- C SETTING UP A SYMMETRICAL COORDINATE SYSTEM. C----------------------------------------------------------------------- DO 10 JG=1,4 SNX(JG)=1.0 10 SNY(JG)=1.0 DO 11 JG=1,6 11 KNK(JG)=1 IF(IFSY.EQ.0)GO TO 20 SNY(3)=-1.0 SNY(4)=-1.0 IS=IS*2 KNK(2)=2 KNK(4)=2 KNK(6)=2 WRITE (ICLA2,1002) 20 IF(IFSX.EQ.0)GO TO 30 SNX(2)=-1.0 SNX(4)=-1.0 IS=IS*2 KNK(2)=2 KNK(3)=2 KNK(4)=2 WRITE (ICLA2,1001) 30 IF(IS.LT.4)GO TO 40 KNK(3)=3 KNK(6)=4 1001 FORMAT(4X,42HTHE FOUNDATION IS SYMMETRICAL ABOUT X-AXIS) 1002 FORMAT(4X,42HTHE FOUNDATION IS SYMMETRICAL ABOUT Y-AXIS) 40 RETURN END C*********************************************************************** SUBROUTINE LUDCMP(A,N,ORDER,TMP) COMPLEX A(N,N),DIV,SUM,TMP(N) INTEGER ORDER(N) DO 10 I=1,N 10 ORDER(I)=I C----FIND PIVOT ELEMENT IN FIRST COLUMN---- CALL APVT(A(1,1),N,ORDER,1) JC=ORDER(1) DIV=A(JC,1) CALL CHKSIZ(DIV) C----FIRST ROW OF U---- DO 20 KCOL=2,N 20 A(JC,KCOL)=A(JC,KCOL)/DIV C----OTHER PARTS OF L AND U---- DO 80 JCOL=2,N-1 JM1=JCOL-1 DO 41 KCOL=1,JM1 KC=ORDER(KCOL) 41 TMP(KCOL)=A(KC,JCOL) C-----DO A COLUMN OF L---- DO 50 IROW=JCOL,N IC=ORDER(IROW) SUM=0.0 DO 40 KCOL=1,JM1 40 SUM=SUM+A(IC,KCOL)*TMP(KCOL) 50 A(IC,JCOL)=A(IC,JCOL)-SUM C-----FIND NEW PIVOT ROW---- CALL APVT(A(1,JCOL),N,ORDER,JCOL) JC=ORDER(JCOL) DIV=A(JC,JCOL) CALL CHKSIZ(DIV) C-----DO A ROW OF U---- DO 70 KCOL=JCOL+1,N SUM=0.0 DO 60 IROW=1,JM1 IC=ORDER(IROW) 60 SUM=SUM+A(JC,IROW)*A(IC,KCOL) 70 A(JC,KCOL)=(A(JC,KCOL)-SUM)/DIV 80 CONTINUE C----LAST ELEMENT OF L---- SUM=0.0 IC=ORDER(N) DO 90 KCOL=1,N-1 KC=ORDER(KCOL) 90 SUM=SUM+A(IC,KCOL)*A(KC,N) A(IC,N)=A(IC,N)-SUM RETURN END C*********************************************************************** SUBROUTINE APVT(A,N,ORDER,JCOL) COMPLEX A(N) INTEGER ORDER(N) IPVT=JCOL JC=ORDER(JCOL) BIG=CABS(A(JC)) DO 10 IROW=JCOL+1,N JC=ORDER(IROW) TEST=CABS(A(JC)) IF(TEST.LT.BIG)GO TO 10 IPVT=IROW BIG=TEST 10 CONTINUE IF(IPVT.EQ.JCOL)RETURN C----SWITCH ROW NUMBERS STORED IN THE ARRAY ORDER---- ISAVE=ORDER(JCOL) ORDER(JCOL)=ORDER(IPVT) ORDER(IPVT)=ISAVE RETURN END C*********************************************************************** SUBROUTINE CHKSIZ(X) COMPLEX X Y=CABS(X) IF(Y.GT.1.E-5)RETURN WRITE (6,100) Y 100 FORMAT(' PIVOT ELEMENT=',E10.3,', IT IS TOO SMALL') RETURN END C*********************************************************************** SUBROUTINE LUSOLV(A,N,B,NI,M,KROWS,ORDER,TMP) COMPLEX A(N,N),B(NI,M),TMP(N),DIV,SUM INTEGER ORDER(N),OFFSET OFFSET=KROWS-1 C----REORDER THE RIGHT HAND SIDE MATRIX---- DO 1 J=1,M DO 2 I=1,N IC=ORDER(I) 2 TMP(I)=B(IC+OFFSET,J) DO 3 I=1,N 3 B(I+OFFSET,J)=TMP(I) 1 CONTINUE C----BACK SUBSTITUTION FOR L---- IC=ORDER(1) DIV=A(IC,1) DO 10 J=1,M 10 B(1+OFFSET,J)=B(1+OFFSET,J)/DIV DO 20 I=2,N IC=ORDER(I) DIV=A(IC,I) DO 21 J=1,M SUM=0.0 DO 22 K=1,I-1 22 SUM=SUM+A(IC,K)*B(K+OFFSET,J) 21 B(I+OFFSET,J)=(B(I+OFFSET,J)-SUM)/DIV 20 CONTINUE C----BACK SUBSTITUTION FOR U---- DO 30 I=N-1,1,-1 IC=ORDER(I) DO 31 J=1,M SUM=0.0 DO 32 K=I+1,N 32 SUM=SUM+A(IC,K)*B(K+OFFSET,J) 31 B(I+OFFSET,J)=B(I+OFFSET,J)-SUM 30 CONTINUE RETURN END C*********************************************************************** SUBROUTINE FOLD(Q,SV,IS,SNX) COMPLEX SV(4),Q(4),TMP1,TMP2 DIMENSION SNX(4) IF(IS.EQ.2)GO TO 100 TMP1=Q(1)+Q(2) TMP2=Q(3)+Q(4) SV(1)=(TMP1+TMP2) SV(3)=(TMP1-TMP2) TMP1=Q(1)-Q(2) TMP2=Q(3)-Q(4) SV(2)=(TMP1+TMP2) SV(4)=(TMP1-TMP2) RETURN 100 SV(1)=(Q(1)+Q(2)) SV(4)=(Q(1)-Q(2)) IF(SNX(4).GT.0.0)GO TO 110 C----SYMMETRY ABOUT Y-AXIS ONLY---- SV(2)=SV(4) SV(3)=SV(1) RETURN C----SYMMETRY ABOUT X-AXIS ONLY---- 110 SV(3)=SV(4) SV(2)=SV(1) RETURN END C*********************************************************************** SUBROUTINE FOLDR(A,NP3F,NTR,NF6,IS,SNX,ISTROW) COMMON/PATTRN/ISW(3,4) COMPLEX A(NTR,NF6),SV(4),Q(4) DIMENSION SNX(4),ILOC(4) INTEGER OFFSET C----FOLDING MATRIX BY ROWS---- NP3=NP3F/IS NP=NP3/3 FACTOR=FLOAT(IS) OFFSET=ISTROW-1 DO 300 ICOMP=1,3 DO 10 KSI=1,IS 10 ILOC(KSI)=ISW(ICOMP,KSI) DO 200 IP=1,NP IR1=(IP-1)*3+ICOMP + OFFSET DO 100 IC=1,NF6 IR=IR1 DO 50 KSI=1,IS Q(KSI)=A(IR,IC) 50 IR=IR+NP3 CALL FOLD(Q,SV,IS,SNX) IR=IR1 DO 60 KSI=1,IS IM=ILOC(KSI) A(IR,IC)=SV(IM)/FACTOR 60 IR=IR+NP3 100 CONTINUE 200 CONTINUE 300 CONTINUE RETURN END C*********************************************************************** SUBROUTINE FOLDL(A,NP3,IS,SNX) COMMON/PATTRN/ISW(3,4) COMPLEX A(NP3,1),SV(4),Q(4) DIMENSION SNX(4),ILOC(4) C----FOLDING MATRIX BY COLUMNS---- NP=NP3/3 NP3F=NP3*IS DO 300 ICOMP=1,3 DO 10 KSI=1,IS 10 ILOC(KSI)=ISW(ICOMP,KSI) DO 200 IP=1,NP IC1=(IP-1)*3+ICOMP DO 100 IR=1,NP3 IC=IC1 DO 50 KSI=1,IS Q(KSI)=A(IR,IC) 50 IC=IC+NP3 CALL FOLD(Q,SV,IS,SNX) IC=IC1 DO 60 KSI=1,IS IM=ILOC(KSI) A(IR,IC)=SV(IM) 60 IC=IC+NP3 100 CONTINUE 200 CONTINUE 300 CONTINUE RETURN END C*********************************************************************** SUBROUTINE UNFOLD(A,NP3F,NTR,NF6,TEMP,IS,SNX,SNY,ISTROW) COMMON/PATTRN/ISW(3,4) COMPLEX A(NTR,NF6),TEMP(NP3F) DIMENSION S(4),SNX(4),SNY(4) INTEGER OFFSET NP3=NP3F/IS NP=NP3/3 OFFSET=ISTROW-1 S(1)=1.0 DO 800 ICOL=1,NF6 DO 100 IPAT=1,NP3F 100 TEMP(IPAT)=0.0 DO 700 KSET=1,IS KSET1=(KSET-1)*NP3 DO 600 ICOMP=1,3 IPAT=ISW(ICOMP,KSET) S(2)=SNX(IPAT) S(3)=SNY(IPAT) S(4)=S(2)*S(3) DO 500 KQUAD=1,IS KQUAD1=(KQUAD-1)*NP3 IPOS=KSET1+ICOMP+OFFSET CALL INPLAC(TEMP(KQUAD1+ICOMP),A(IPOS,ICOL),NP3,S(KQUAD)) 500 CONTINUE 600 CONTINUE 700 CONTINUE DO 200 IPAT=1,NP3F 200 A(IPAT+OFFSET,ICOL)=TEMP(IPAT) 800 CONTINUE RETURN END C*********************************************************************** SUBROUTINE INPLAC(TEMP,A,NP3,SKQ) COMPLEX TEMP(1),A(1) IF(SKQ.LT.0.0)GO TO 200 DO 100 I=1,NP3,3 100 TEMP(I)=TEMP(I)+A(I) RETURN 200 DO 300 I=1,NP3,3 300 TEMP(I)=TEMP(I)-A(I) RETURN END C***********************************************************************