#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'subs.f' <<'END_OF_FILE' X PROGRAM MAIN X C THIS PROGRAM WILL TEST THE ABILITY OF A SOFTWARE/HARDWARE C SYSTEM IN REGARDS TO PARALLEL CONSTRUCTS. C THE PROGRAM CONSISTS OF THE MAIN PROGRAM, DRIVERS THAT C CALL AND TIME THE SUBROUTINE TO BE TESTED AND A COLLECTION C OF SUBROUTINES FROM VARIOUS PROGRAMS. THESE SUBROUTINES HAVE BEEN C ASSEMBLED FROM A VARIETY OF PROGRAMS AND AS SUCH THEY HAVE BEEN C TAKEN OUT OF THEIR PROPER CONTEXT. THEY HAVE BEEN MODIFIED SO C THAT THEY WILL EXECUTE IN A STAND ALONE MODE. C THE PARAMETER 'NOROUT' SETS UP THE NUMBER OF SUBROUTINES C TO TEST. IF ROUTINES ARE ADDED THIS PARAMETER MUST BE INCREASED C TO THE CORRECT NUMBER. THE 'WTIMES' ARRAY CONTAINS THE WALL CLOCK C TIMES FOR EACH SUBROUTINE THAT IS TESTED. C THE 'IANS' ARRAY CONTAINS A O IF THE ROUTINE CALCULATED TO A C PREDETERMINED RESULT AND A 1 IF IT DID NOT OBTAIN THE DESIRED RESULT. C 'ITERNO IS A LOOP COUNTER THAT ALLOWS THE AMOUNT OF CPU TIME C GENERATED TO BE VARIED. EACH DRIVER ROUTINE GENERATES AROUND .3 OF A C CPU SECOND ON A CRAY XMP. THE DRIVER ALSO RETURNS THE CPU AND WALL C CLOCK TIME FOR EACH ROUTINE. THERE IS A DRIVER FOR EACH ROUTINE C THAT IS CALLED. THE 'IEX' ARRAY IS A SET OF FLAGS WHICH STATES C WHAT ROUTINES TO TEST. A 1 IN IEX(1) WILL CALL DRIVER 1. A 0 WILL C SKIP CALLING THAT ROUTINE. INSIDE OF EACH DRIVER THE ROUTINE TO BE C TESTED MAY BE CALLED MORE THAN ONCE IN ORDER TO GENERATE ENOUGH C CPU TIME. VARIABLE 'PRTALL' CONTROLS THE PRINTING OF INDIVIDUAL C TIMINGS FOR EACH ROUTINE. IF IT IS 'TRUE' , THEN THE TIMES FOR C EACH ROUTINE WILL BE PRINTED. IF 'FALSE' ONLY THE TOTAL TIMES WILL C BE PRINTED. VARIABLE 'ITERNO' CONTROLS THE NUMBER OF TIMES THAT C THE PROGRAMS WILL BE EXECUTED IN AN OUTER LOOP. THIS ALLOWS THE C GENERATION OF MORE OR LESS TIME. THE MFLOPS ARE CALCULATED BY C DIVIDING THE PREDETERMINED FLOPS BY THE WALL CLOCK TIME. X X PARAMETER ( NOROUT = 36 , NWORDS = 2100000) X X COMMON /DATA/DUMMY(NWORDS) X COMMON /ERROR/ RELERR X COMMON/WT/ WTIMES(NOROUT) X COMMON /RESULTS/ IANS(NOROUT) X X DIMENSION IEX(NOROUT),FLOPS(NOROUT) X X INTEGER FLOPS,FLOPST X REAL MFLOPS,MFLOPA X LOGICAL PRTALL X X DATA PRTALL/.TRUE./ X DATA IEX/NOROUT*1/ X DATA ITERNO/1/ X DATA FLOPS/ X * 7303347, 5660969, 126307648, 270513504, 17668534, 173757577, X * 291647851, 24518837, 11256093, 10452045, 28981773, 101543326, X * 10304303, 221141203, 12815237, 33089704, 243753086, 10121093, X * 190660062, 6913440, 23763905, 130846623, 4904591, 196496290, X * 22563098, 12913459, 1323307, 13717990, 180614247, 334903167, X * 136813658, 122346352, 2528014, 27661012, 36070850, 16585440/ X DATA CONST/1.E-06/ X X RELERR = 1.E-9 X DO 1000 I = 1, NOROUT X WTIMES(I) = 0 X IANS(I) = 0 X 1000 CONTINUE X X DO 1200 I = 1, ITERNO X IF (IEX(1) .NE. 0) CALL D1 X IF (IEX(2) .NE. 0) CALL D2 X IF (IEX(3) .NE. 0) CALL D3 X IF (IEX(4) .NE. 0) CALL D4 X IF (IEX(5) .NE. 0) CALL D5 X IF (IEX(6) .NE. 0) CALL D6 X IF (IEX(7) .NE. 0) CALL D7 X IF (IEX(8) .NE. 0) CALL D8 X IF (IEX(9) .NE. 0) CALL D9 X IF (IEX(10) .NE. 0) CALL D10 X IF (IEX(11) .NE. 0) CALL D11 X IF (IEX(12) .NE. 0) CALL D12 X IF (IEX(13) .NE. 0) CALL D13 X IF (IEX(14) .NE. 0) CALL D14 X IF (IEX(15) .NE. 0) CALL D15 X IF (IEX(16) .NE. 0) CALL D16 X IF (IEX(17) .NE. 0) CALL D17 X IF (IEX(18) .NE. 0) CALL D18 X IF (IEX(19) .NE. 0) CALL D19 X IF (IEX(20) .NE. 0) CALL D20 X IF (IEX(21) .NE. 0) CALL D21 X IF (IEX(22) .NE. 0) CALL D22 X IF (IEX(23) .NE. 0) CALL D23 X IF (IEX(24) .NE. 0) CALL D24 X IF (IEX(25) .NE. 0) CALL D25 X IF (IEX(26) .NE. 0) CALL D26 X IF (IEX(27) .NE. 0) CALL D27 X IF (IEX(28) .NE. 0) CALL D28 X IF (IEX(29) .NE. 0) CALL D29 X IF (IEX(30) .NE. 0) CALL D30 X IF (IEX(31) .NE. 0) CALL D31 X IF (IEX(32) .NE. 0) CALL D32 X IF (IEX(33) .NE. 0) CALL D33 X IF (IEX(34) .NE. 0) CALL D34 X IF (IEX(35) .NE. 0) CALL D35 X IF (IEX(36) .NE. 0) CALL D36 X 1200 CONTINUE X X IFAILS = 0 X DO 1600 IDNUM = 1, NOROUT X IF (IEX(IDNUM) .EQ. 0) GO TO 1400 X IF (IANS(IDNUM) .NE. 0) IFAILS = IFAILS + 1 X 1400 CONTINUE X 1600 CONTINUE X IF (IFAILS .EQ. 0) THEN X PRINT 12 X ELSE X PRINT 14 X DO 2000 IDNUM = 1, NOROUT X IF (IEX(IDNUM) .EQ. 0) GO TO 1800 X IF (IANS(IDNUM) .NE. 0) PRINT 10, IDNUM X 1800 CONTINUE X 2000 CONTINUE X ENDIF X X WLLTME = 0 X FLOPST = 0 X EXECS = 0 X GM = 1. X DO 2400 IDNUM = 1, NOROUT X IF (IEX(IDNUM) .EQ. 0) GO TO 2200 X WLLTME = WLLTME + WTIMES(IDNUM) X IF (ITERNO .NE. 1) FLOPS(IDNUM) = FLOPS(IDNUM)*ITERNO X FLOPST = FLOPST + FLOPS(IDNUM) X MFLOPS = (FLOPS(IDNUM)/WTIMES(IDNUM))*CONST X IF (PRTALL) PRINT 16, IDNUM, WTIMES(IDNUM), MFLOPS X GM = GM * MFLOPS X EXECS = EXECS + 1. X 2200 CONTINUE X 2400 CONTINUE X X MFLOPA = (FLOPST/WLLTME)*CONST X X GM = GM**(1./EXECS) X PRINT 18, WLLTME, MFLOPA X PRINT 20, GM X STOP X X X 10 FORMAT(1H ,'RESULTS FROM SUBROUTINE',I3,' WERE INVALID') X 12 FORMAT(1H ,'ALL RESULTS WERE VALID') X 14 FORMAT(1H ,'NOT ALL RESULTS WERE VALID. ') X 16 FORMAT(1H ,'SUB',I2,2X,' WALL CLOCK SECONDS ',F6.2,2X X 1 ,'MFLOPS ',F7.1) X 18 FORMAT(1H ,/,' TOTAL WALL CLOCK SECONDS ',F8.3,/, X 1 ' AVERAGE MFLOPS ',F8.3) X 20 FORMAT(1H ,/,' GEOMETRIC MEAN ',F8.3) X END X X SUBROUTINE CALLP0(ID,SUB) X PARAMETER (CLOCK=6.0E-09) X COMMON /OVERH/ OVER1,OVER2,OVER3 X COMMON/WT/ WTIMES(1) X INTEGER WALLS, WALLE X WALLS = RTC() X CALL SUB X WALLE = RTC() X GO TO 1000 X X ENTRY CALLP1 (ID, SUB, A) X WALLS = RTC() X CALL SUB (A) X WALLE = RTC() X GO TO 1000 X X ENTRY CALLP2 (ID, SUB, A, B) X WALLS = RTC() X CALL SUB (A, B) X WALLE = RTC() X GO TO 1000 X X ENTRY CALLP3 (ID, SUB, A, B, C) X WALLS = RTC() X CALL SUB (A, B, C) X WALLE = RTC() X GO TO 1000 X X ENTRY CALLP4 (ID, SUB, A, B, C, D) X WALLS = RTC() X CALL SUB (A, B, C, D) X WALLE = RTC() X GO TO 1000 X X ENTRY CALLP5 (ID, SUB, A, B, C, D, E) X WALLS = RTC() X CALL SUB (A, B, C, D, E) X WALLE = RTC() X GO TO 1000 X X ENTRY CALLP6 (ID, SUB, A, B, C, D, E, F) X WALLS = RTC() X CALL SUB (A, B, C, D, E, F) X WALLE = RTC() X GO TO 1000 X X ENTRY CALLP7 (ID, SUB, A, B, C, D, E, F, G) X WALLS = RTC() X CALL SUB (A, B, C, D, E, F, G) X WALLE = RTC() X GO TO 1000 X X ENTRY CALLP8 (ID, SUB, A, B, C, D, E, F, G, H) X WALLS = RTC() X CALL SUB (A, B, C, D, E, F, G, H) X WALLE = RTC() X X 1000 CONTINUE X WTIMES(ID) = WTIMES(ID) + (FLOAT((WALLE-WALLS))*CLOCK) X RETURN X END X X SUBROUTINE INIT(VAR,CON,N,INC) X COMPLEX CVAR(N),CCON X DIMENSION VAR(N),IVAR(N) X DO 1000 I = 1, N, INC X VAR(I) = CON X 1000 CONTINUE X RETURN X X ENTRY IINIT (IVAR, ICON, N, INC) X DO 1200 I = 1, N, INC X IVAR(I) = ICON X 1200 CONTINUE X RETURN X X ENTRY CINIT (CVAR, CCON, N, INC) X DO 1400 I = 1, N, INC X CVAR(I) = CCON X 1400 CONTINUE X RETURN X END X X SUBROUTINE RSUM(ID,VAR,N,RESULT) X COMPLEX CVAR(N),CRSULT,CRES X DIMENSION VAR(N),IVAR(N) X COMMON /ERROR/ RELERR X COMMON /RESULTS/ IANS(1) X X SUM = 0 X IANSWR = 0 X DO 1000 I = 1, N X SUM = SUM + VAR(I) X 1000 CONTINUE X DIFF = ABS(RESULT-SUM) X ESUM = ABS(RELERR*SUM) X IF (DIFF .GT. ESUM) IANSWR = 1 X IF (IANSWR .NE. 0) PRINT 10, SUM, RESULT X 10 FORMAT(' CALC ',E20.14,/,' EXPC ',E20.14) X GO TO 1600 X X ENTRY IISUM (ID, IVAR, N, IRSULT) X IRES = 0 X IANSWR = 0 X DO 1200 I = 1, N X IRES = IRES + IVAR(I) X 1200 CONTINUE X DIFF = ABS(IRSULT-IRES) X ESUM = ABS(RELERR*IRES) X IF (DIFF .GT. ESUM) IANSWR = 1 X GO TO 1600 X X ENTRY CCSUM (ID, CVAR, N, CRSULT) X CRES = 0 X IANSWR = 0 X DO 1400 I = 1, N X CRES = CRES + CVAR(I) X 1400 CONTINUE X DIFF = ABS(CRSULT-CRES) X ESUM = ABS(RELERR*CRES) X IF (DIFF .GT. ESUM) IANSWR = 1 X 1600 CONTINUE X IANS(ID) = IANS(ID) + IANSWR X RETURN X END X SUBROUTINE D1 X PARAMETER (M1=400000,N=1000,M=100) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB001 X DIMENSION CON(4) X DATA ANS /.23465302666590E+3/ X DATA CON(1)/.0014/,CON(2)/.005678/,CON(3)/.0010/ X DATA CON(4)/.001234/ X DO 1200 I = 1, 12 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP2 (1, SUB001, M, N) X 1200 CONTINUE X CALL RSUM (1, DUM1, M*N, ANS) X RETURN X END X SUBROUTINE D2 X PARAMETER (N=10,M=10) X PARAMETER (M1=N*M) X COMMON /DATA/ DUM1(M),DUM2(M1),DUM3(M1),DUM4(M1),DUM5(M1) X EXTERNAL SUB002 X DIMENSION CON(4),A2(N,M) X DATA ANS /.35722539999999E+2/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 5000 X DO 1000 J = 1, 4 X CALL INIT (DUM5(J),CON(J),N*M,4) X CALL INIT (DUM4(J),CON(J),N*M,4) X CALL INIT (DUM3(J),CON(J),N*M,4) X CALL INIT (DUM2(J),CON(J),N*M,4) X CALL INIT (DUM1(J),CON(J),N*M,4) X 1000 CONTINUE X CALL CALLP7 (2, SUB002, M, N, DUM1, DUM2, DUM3, DUM4, DUM5) X 1200 CONTINUE X CALL RSUM (2, DUM2, M*N, ANS) X RETURN X END X SUBROUTINE D3 X PARAMETER (N=100,M=100) X PARAMETER (M1=M*N) X COMMON /DATA/ DUM1(M),DUM2(M1),DUM3(M1),DUM4(M1),DUM5(M1) X EXTERNAL SUB003 X DIMENSION CON(4) X DATA ANS /.35722539999987E+4/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 1400 X DO 1000 J = 1, 4 X CALL INIT (DUM5(J),CON(J),N*M,4) X CALL INIT (DUM4(J),CON(J),N*M,4) X CALL INIT (DUM3(J),CON(J),N*M,4) X CALL INIT (DUM2(J),CON(J),N*M,4) X CALL INIT (DUM1(J),CON(J),N*M,4) X 1000 CONTINUE X CALL CALLP7 (3, SUB003, M, N, DUM1, DUM2, DUM3, DUM4, DUM5) X 1200 CONTINUE X CALL RSUM (3, DUM2, M*N, ANS) X RETURN X END X SUBROUTINE D4 X PARAMETER (N=500,M=1000) X PARAMETER (M1=M*N) X COMMON /DATA/ DUM1(M),DUM2(M1),DUM3(M1),DUM4(M1),DUM5(M1) X EXTERNAL SUB004 X DIMENSION CON(4) X DATA ANS /.17861269999898E+6/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 60 X DO 1000 J = 1, 4 X CALL INIT (DUM5(J),CON(J),M1,4) X CALL INIT (DUM4(J),CON(J),M1,4) X CALL INIT (DUM3(J),CON(J),M1,4) X CALL INIT (DUM2(J),CON(J),M1,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP7 (4, SUB004, M, N, DUM1, DUM2, DUM3, DUM4, DUM5) X 1200 CONTINUE X CALL RSUM (4, DUM2, M*N, ANS) X RETURN X END X SUBROUTINE D5 X PARAMETER (M1=21*101*101) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB005 X DIMENSION CON(4) X DATA ANS /-.42552817357591E-3/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 20 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP1 (5, SUB005, SUM) X 1200 CONTINUE X CALL RSUM (5, SUM, 1, ANS) X RETURN X END X SUBROUTINE D6 X PARAMETER (M1=6*100*100*10+9*10*100+5*100) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB006 X DIMENSION CON(4) X DATA ANS /.80019108524855E+2/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 42 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP0 (6, SUB006) X 1200 CONTINUE X CALL RSUM (6, DUM1, 200, ANS) X RETURN X END X SUBROUTINE D7 X PARAMETER (M1=22*101*101) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB007 X DIMENSION CON(4) X DATA ANS /.18733898409327E+5/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 350 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP0 (7, SUB007) X 1200 CONTINUE X CALL RSUM (7, DUM1, 3*101*101, ANS) X RETURN X END X SUBROUTINE D8 X PARAMETER (M1=100*1000,M2=1000) X COMMON /DATA/ DUM1(M1),DUM2(M2),DUM3(M2) X EXTERNAL SUB008 X DIMENSION CON(4) X DATA ANS /.18724499999957E+5/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 8 X DO 1000 J = 1, 4 X CALL INIT (DUM3(J),CON(J),M2,4) X CALL INIT (DUM2(J),CON(J),M2,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP6 (8, SUB008, 1000, 100, DUM2, DUM3, DUM1, IDUM) X 1200 CONTINUE X CALL RSUM (8, DUM1, M1, ANS) X RETURN X END X SUBROUTINE D9 X PARAMETER (M1=27000+6*30,M2=27000*3) X COMMON /DATA/ DUM1(M1),DUM2(M2) X EXTERNAL SUB009 X DIMENSION CON(4) X DATA ANS /.79003609513564E+3/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 100 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL INIT (DUM2, 0, M2, 1) X CALL CALLP4 (9, SUB009, .5, 1.5, 3.14, DUM2) X 1200 CONTINUE X CALL RSUM (9, DUM2, M2, ANS) X RETURN X END X SUBROUTINE D10 X PARAMETER (M1=1600000,M2=12*12*24*6*24) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB010 X DIMENSION CON(4),IDUM(4) X DATA IDUM/4*1/ X DATA ANS /-.10353899519837E+6/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 2 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP3 (10, SUB010, 1, 1, IDUM) X 1200 CONTINUE X CALL RSUM (10, DUM1, M2, ANS) X RETURN X END X SUBROUTINE D11 X PARAMETER (M1=1000+1000*1000) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB011 X DIMENSION CON(4) X DATA ANS /.11420219304919E+0/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 7 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP0 (11, SUB011) X 1200 CONTINUE X CALL RSUM (11, DUM1, 1000, ANS) X RETURN X END X SUBROUTINE D12 X PARAMETER (II=101,JJ=101,KK=26) X PARAMETER (IJK=II*JJ*KK) X PARAMETER (M1=IJK+260) X COMMON /DATA/ DUM1(M1),DUM2(IJK) X EXTERNAL SUB012 X DIMENSION CON(4) X DATA ANS /.37440876253074E+3/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 7 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL INIT (DUM2, .00014576, IJK, 1) X CALL CALLP1 (12, SUB012, DUM2) X 1200 CONTINUE X CALL RSUM (12, DUM2, IJK, ANS) X RETURN X END X SUBROUTINE D13 X PARAMETER (M1=4*1000*100+5*1000+10*100) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB013 X DIMENSION CON(4) X DATA ANS /.18872364107147E+1/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 6 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP0 (13, SUB013) X 1200 CONTINUE X CALL RSUM (13, DUM1, 100, ANS) X RETURN X END X SUBROUTINE D14 X PARAMETER (M1=1900000) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB014 X DIMENSION CON(4) X DATA ANS /.47934719999997E+3/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 3 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP0 (14, SUB014) X 1200 CONTINUE X CALL RSUM (14, DUM1, 48*72*32, ANS) X RETURN X END X SUBROUTINE D15 X PARAMETER (M1=110000,M=10000,NM=172) X COMMON /DATA/ DUM1(M1),DUM2(M),IDUM3(NM),IDUM4(NM) X EXTERNAL SUB015 X DIMENSION CON(4) X DATA ANS /.75126427075617E-1/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/3.1234/ X DO 1200 I = 1, 2 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL INIT (DUM2, 0, M, 1) X CALL IINIT (IDUM3, 0, NM, 1) X CALL IINIT (IDUM4, 0, NM, 1) X CALL CALLP3 (15, SUB015, DUM2, IDUM3, IDUM4) X 1200 CONTINUE X CALL RSUM (15, DUM2, M, ANS) X RETURN X END X SUBROUTINE D16 X COMPLEX DUM1,DUM2,ANS,CON X PARAMETER (M1=10000) X COMMON /DATA/ I1,I2,I3,I4,I5,I6,I7,DUM1(M1),DUM2(2) X EXTERNAL SUB016 X DIMENSION CON(4) X DATA ANS /(-.36348926727874E-14,-.24624692117746E-14)/ X DATA CON(1)/(.14,.14)/,CON(2)/(.5678,.5678)/ X * ,CON(3)/(.001,.001)/,CON(4)/(.1234,.1234)/ X I1 = 5000 X I2 = 1 X I3 = 200 X I4 = 1 X I5 = 0 X I6 = 1 X I7 = 1 X DO 1000 J = 1, 4 X CALL CINIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X DUM2(1) = 1. X DUM2(2) = 1. X 1200 CONTINUE X CALL CALLP2 (16, SUB016, DUM1, DUM2) X CALL CCSUM (16, DUM1, 5000, ANS) X RETURN X END X SUBROUTINE D17 X PARAMETER (NNX=44) X PARAMETER (NNY=44) X PARAMETER (NNZ=4) X PARAMETER (MXC=11) X PARAMETER (MXP=3) X PARAMETER (MXW=5) X PARAMETER (NXYN=NNX*NNY) X PARAMETER (NBL=NXYN*NNZ) X PARAMETER (NBLW=NBL+MXW) X PARAMETER (M1=1589000,M2=60) X COMMON /DATA/ DUM1(M1) X COMMON /INTS/IDUM1(M2) X EXTERNAL SUB017 X DIMENSION CON(4) X DATA ANS /.29400934903364E+11/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 4 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL IINIT (IDUM1, 1, M2, 1) X CALL CALLP0 (17, SUB017) X 1200 CONTINUE X CALL RSUM (17, DUM1, NBLW*MXC, ANS) X RETURN X END X SUBROUTINE D18 X PARAMETER (M1=700000,M2=100000) X COMMON /DATA/ DUM1(M1),DUM2(M2),DUM3(M2),DUM4(M2) X EXTERNAL SUB018 X DIMENSION CON(4) X DATA ANS /.25530009809248E+5/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.19/,CON(4)/.1234/ X DO 1200 I = 1, 100 X DO 1000 J = 1, 4 X CALL INIT (DUM4(J),CON(J),M2,4) X CALL INIT (DUM3(J),CON(J),M2,4) X CALL INIT (DUM2(J),CON(J),M2,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP4 (18, SUB018, DUM1, DUM2, DUM3, DUM4) X 1200 CONTINUE X CALL RSUM (18, DUM4, M2, ANS) X RETURN X END X SUBROUTINE D19 X PARAMETER (M1=688000) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB019 X DIMENSION CON(4) X DATA ANS /.58082299171138E+5/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 70 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP0 (19, SUB019) X 1200 CONTINUE X CALL RSUM (19, DUM1, M1, ANS) X RETURN X END X SUBROUTINE D20 X PARAMETER (M1=29791*9+100*5,M2=29791*5,M3=100) X COMMON /DATA/ DUM1(M1),DUM2(M2),IDUM3(M3),IDUM4(M3) X EXTERNAL SUB020 X DIMENSION CON(4) X DATA ANS /.46771869999876E+6/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL INIT (DUM2, 3.14, M2, 1) X CALL IINIT (IDUM3, 1, M3, 1) X CALL IINIT (IDUM4, 1, M3, 1) X 1200 CONTINUE X CALL CALLP3 (20, SUB020, DUM2, IDUM3, IDUM4) X CALL RSUM (20, DUM2, M2, ANS) X RETURN X END X SUBROUTINE D21 X PARAMETER (M1=100000) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB021 X DIMENSION CON(4) X DATA ANS /.35916589036686E+11/ X DATA CON(1)/.4/,CON(2)/.678/,CON(3)/.9/, X * CON(4)/.234/ X DO 1200 I = 1, 3 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP4 (21, SUB021, 100, 100, 100, 0) X 1200 CONTINUE X CALL RSUM (21, DUM1, 6000, ANS) X RETURN X END X SUBROUTINE D22 X PARAMETER (M1=625175,M2=31250) X COMMON /DATA/ DUM1(M1),DUM2(M2) X EXTERNAL SUB022 X DIMENSION CON(4) X DATA ANS /.30710253481606E+8/ X DATA CON(1)/1.4/,CON(2)/5.678/,CON(3)/1.9/,CON(4)/1.234/ X DO 1200 I = 1, 14 X DO 1000 J = 1, 4 X CALL INIT (DUM2(J),CON(J),M2,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP4 (22, SUB022, DUM2, 1, 150, 5) X 1200 CONTINUE X CALL RSUM (22, DUM2, M2, ANS) X RETURN X END X SUBROUTINE D23 X PARAMETER (M1=513000) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB023 X DIMENSION CON(4) X DATA ANS /-.34586383403291E-2/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 1 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP1 (23, SUB023, SUM) X 1200 CONTINUE X CALL RSUM (23, SUM, 1, ANS) X RETURN X END X SUBROUTINE D24 X PARAMETER (M1=10200,M2=10000) X COMMON /DATA/ DUM1(M1),DUM2(M2),DUM3(M2),DUM4(M2) X EXTERNAL SUB024 X DIMENSION CON(4) X DATA ANS /-.20207560939631E+3/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.19/,CON(4)/.1234/ X DO 1200 I = 1, 21 X DO 1000 J = 1, 4 X CALL INIT (DUM2(J),CON(J),M2,4) X CALL INIT (DUM3(J),CON(J),M2,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL INIT (DUM4, .0001, M2, 1) X CALL CALLP3 (24, SUB024, DUM2, DUM3, DUM4) X 1200 CONTINUE X CALL RSUM (24, DUM2, M2*2, ANS) X RETURN X END X SUBROUTINE D25 X PARAMETER (M1=50000,IJ=28) X COMMON /DATA/ DUM1(M1),DUM2(IJ) X EXTERNAL SUB025 X DIMENSION CON(4) X DATA ANS /.26728212215467E+6/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 3 X DO 1000 J = 1, 4 X CALL INIT (DUM2(J),CON(J),IJ,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP1 (25, SUB025, DUM2) X 1200 CONTINUE X CALL RSUM (25, DUM2, IJ, ANS) X RETURN X END X SUBROUTINE D26 X PARAMETER (M1=36,M2=47000,M3=44) X COMPLEX DUM1,CCON X COMMON /DATA/ DUM1(M1),DUM2(M2),DUM3(M3) X EXTERNAL SUB026 X DIMENSION CON(4) X DATA ANS /-.78506753418025E+2/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.19/,CON(4)/.1234/ X DO 1200 I = 1, 12 X DO 1000 J = 1, 4 X CALL INIT (DUM2(J),CON(J),M2,4) X CCON = CON(J) X CALL CINIT (DUM1(J),CCON,M1,4) X 1000 CONTINUE X CALL INIT (DUM3, 0, M3, 1) X CALL CALLP1 (26, SUB026, DUM3) X 1200 CONTINUE X CALL RSUM (26, DUM3, M3, ANS) X RETURN X END X SUBROUTINE D27 X PARAMETER (M1=100*100*100,IJ=40) X COMMON /DATA/ DUM1(M1),DUM2(IJ) X EXTERNAL SUB027 X DIMENSION CON(4) X DATA ANS /.21304319999997E+5/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 16 X DO 1000 J = 1, 4 X CALL INIT (DUM2(J),CON(J),IJ,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP0 (27, SUB027) X 1200 CONTINUE X CALL RSUM (27, DUM1, M1, ANS) X RETURN X END X SUBROUTINE D28 X PARAMETER (M1=200500) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB028 X DIMENSION CON(4) X DATA ANS /.20164390072275E+5/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 2 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP0 (28, SUB028) X 1200 CONTINUE X CALL RSUM (28, DUM1, 1000*100, ANS) X RETURN X END X SUBROUTINE D29 X PARAMETER (M1=44000) X COMMON /DATA/ DUM1(M1),DUM2(25) X EXTERNAL SUB029 X DIMENSION CON(4) X DATA ANS /.14228764570541E+6/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X Do 1200 I = 1, 5 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL INIT(DUM2,0.0,25,1) X CALL CALLP1 (29, SUB029, DUM2) X 1200 CONTINUE X CALL RSUM (29, DUM2, 25, ANS) X RETURN X END X SUBROUTINE D30 X PARAMETER (M1=800000,I2=1000,J2=100) X COMMON /DATA/ DUM1(M1),DUM2(J2) X EXTERNAL SUB030 X DIMENSION CON(4) X DATA ANS /.66713894216008E+5/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 200 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP3 (30, SUB030, J2, I2, DUM2) X 1200 CONTINUE X CALL RSUM (30, DUM2, J2, ANS) X RETURN X END X SUBROUTINE D31 X PARAMETER (LS=1,LE=101,KS=1,KE=129,KD=129,LD=101,MD=129) X PARAMETER (M1=MD*25,M2=LD*KD*125,M3=LD*75,M4=LD*5,M5=LD*KD*5) X PARAMETER (M6=M2+M3+M4+M5) X COMMON /DATA/ DUM1(M1),DUM2(M6) X EXTERNAL SUB031 X DATA ANS / .32234316520824E+02 / X DATA CON1 /.01 /, CON2 /1.0/ X DO 1000 I = 1, 6 X CALL INIT (DUM1, CON2, M1, 1) X CALL INIT (DUM2, CON1, M6, 1) X CALL CALLP0 (31, SUB031) X 1000 CONTINUE X CALL RSUM (31, DUM2, M5, ANS) X RETURN X END X SUBROUTINE D32 X PARAMETER (JMAX=120,KMAX=23,LMAX=30) X PARAMETER (M2=JMAX*KMAX*LMAX*17,M3=JMAX*JMAX*4) X PARAMETER (M1=JMAX*KMAX*LMAX*6,M4=M2+M3+5) X COMMON /DATA/ DUM1(M1),DUM2(M4) X EXTERNAL SUB032 X DATA ANS /.12279940726070E+6/ X DATA CON1/3.14/,CON2/.5678/ X DO 1000 I = 1, 5 X CALL INIT (DUM1, CON2, M1, 1) X CALL INIT (DUM2, CON1, M4, 1) X CALL CALLP0 (32, SUB032) X 1000 CONTINUE X CALL RSUM (32, DUM1, M1, ANS) X RETURN X END X SUBROUTINE D33 X PARAMETER (M1=15000,I2=1000) X COMMON /DATA/ DUM1(M1),JM1(4*I2),JM2(I2),JM3(I2),JM4(I2) X COMMON /DATA/ JM5(I2),JM6(I2) X EXTERNAL SUB033 X DIMENSION CON(4) X DATA ANS /.36315788655113E+4/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 15 X DO 1000 J = 1, 4 X CALL IINIT (JM1(J),J,9000,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP7 (33, SUB033, JM1, JM2, JM3, JM4, JM5, JM6, I2) X 1200 CONTINUE X CALL RSUM (33, DUM1, I2*5, ANS) X RETURN X END X SUBROUTINE D34 X PARAMETER (M1=37*100*100) X COMMON /DATA/ DUM1(M1) X EXTERNAL SUB034 X DIMENSION CON(4) X DATA ANS /.61065011581289E+6/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 10 X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP2 (34, SUB034, 100, 100) X 1200 CONTINUE X CALL RSUM (34, DUM1, 100*100, ANS) X RETURN X END X SUBROUTINE D35 X PARAMETER (M1=500000,I2=96,J2=24,K2=24) X COMMON /DATA/ DUM1(M1),DUM2(M1),DUM3(M1),DUM4(M1) X EXTERNAL SUB035 X DIMENSION CON(4) X DATA ANS /.13346223580218E+4/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1200 I = 1, 2 X DO 1000 J = 1, 4 X CALL INIT (DUM4(J),CON(J),M1,4) X CALL INIT (DUM3(J),CON(J),M1,4) X CALL INIT (DUM2(J),CON(J),M1,4) X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP4 (35, SUB035, I2, J2, K2, RTRMS) X 1200 CONTINUE X CALL RSUM (35, RTRMS, 1, ANS) X RETURN X END X SUBROUTINE D36 X PARAMETER (M1=241000,N=20,M=80) X COMMON /DATA/ DUM1(M1),DUM2(N*M) X EXTERNAL SUB036 X DIMENSION CON(4) X DATA ANS /.32343671293395E+4/ X DATA CON(1)/.14/,CON(2)/.5678/,CON(3)/.0010/,CON(4)/.1234/ X DO 1000 J = 1, 4 X CALL INIT (DUM1(J),CON(J),M1,4) X 1000 CONTINUE X CALL CALLP3 (36, SUB036, N, M, DUM2) X CALL RSUM (36, DUM2, N*M, ANS) X RETURN X END X SUBROUTINE SUB001(M,N) X PARAMETER (MM=100,NN=1000) X COMMON /DATA/ A(MM,NN) X COMMON /DATA/ X(MM,NN),Y(MM,NN),Z(MM,NN) X X S = .0001 X DO 1400 I = 1, M X DO 1000 J = 1, N X A(I,J) = A(I,J) + 1 X IF (A(I,J) .EQ. 0) GO TO 1200 X 1000 CONTINUE X 1200 CONTINUE X 1400 CONTINUE X X DO 1800 I = 1, M X DO 1600 J = 1, N X A(I,J) = 1 X 1600 CONTINUE X A(I,I) = A(I,I) + 1 X 1800 CONTINUE X X DO 2200 I = 1, M X DO 2000 J = 1, N X A(I,J) = S X Y(I,J) = (X(I,J)+Z(I,J))*S X S = Y(I,J) + X(I,J) X 2000 CONTINUE X 2200 CONTINUE X X DO 2600 I = 1, M X DO 2400 J = 1, N X X(I,J) = Y(I,J) X Y(I,J) = X(I,J)*Z(I,J) X 2400 CONTINUE X 2600 CONTINUE X X RETURN X END X SUBROUTINE SUB002(M,N,A1,A2,B2,C2,D2) X DIMENSION A2(M,N) X REAL A1(M) X REAL B2(M,N) X REAL C2(M,N) X REAL D2(M,N) X DATA S1/0/ C C DOUBLY DIMENSIONED ARRAY - SIMPLE ASSIGNMENT C X DO 1200 I = 1, N X DO 1000 J = 1, M X A2(J,I) = B2(J,I) + C2(J,I) X 1000 CONTINUE X 1200 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH SCALAR MULTIPLIER C X DO 1600 I = 1, N X DO 1400 J = 1, M X A2(J,I) = B2(J,I) + S1*C2(J,I) X 1400 CONTINUE X 1600 CONTINUE C C DOUBLY DIMENSIONED ARRAY - WITH REUSED VALUE C X DO 2000 I = 1, N X DO 1800 J = 1, M X A2(J,I) = B2(J,I) + C2(J,I) X B2(J,I) = A2(J,I) + D2(J,I) X 1800 CONTINUE X 2000 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH PRIVATE SCALAR TEMPORARY C X DO 2400 I = 1, N X DO 2200 J = 1, M X A06 = B2(J,I) + C2(J,I) X A2(J,I) = A06*D2(J,I) X 2200 CONTINUE X 2400 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH PRIVATE ARRAY TEMPORARY C X DO 2800 I = 1, N X DO 2600 J = 1, M X A1(J) = B2(J,I) + C2(J,I) X A2(J,I) = A1(J)*D2(J,I) X 2600 CONTINUE X 2800 CONTINUE X RETURN X END X SUBROUTINE SUB003(M,N,A1,A2,B2,C2,D2) X DIMENSION A2(M,N) X REAL A1(M) X REAL B2(M,N) X REAL C2(M,N) X REAL D2(M,N) X DATA S1/0/ C C DOUBLY DIMENSIONED ARRAY - SIMPLE ASSIGNMENT C X DO 1200 I = 1, N X DO 1000 J = 1, M X A2(J,I) = B2(J,I) + C2(J,I) X 1000 CONTINUE X 1200 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH SCALAR MULTIPLIER C X DO 1600 I = 1, N X DO 1400 J = 1, M X A2(J,I) = B2(J,I) + S1*C2(J,I) X 1400 CONTINUE X 1600 CONTINUE C C DOUBLY DIMENSIONED ARRAY - WITH REUSED VALUE C X DO 2000 I = 1, N X DO 1800 J = 1, M X A2(J,I) = B2(J,I) + C2(J,I) X B2(J,I) = A2(J,I) + D2(J,I) X 1800 CONTINUE X 2000 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH PRIVATE SCALAR TEMPORARY C X DO 2400 I = 1, N X DO 2200 J = 1, M X A06 = B2(J,I) + C2(J,I) X A2(J,I) = A06*D2(J,I) X 2200 CONTINUE X 2400 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH PRIVATE ARRAY TEMPORARY C X DO 2800 I = 1, N X DO 2600 J = 1, M X A1(J) = B2(J,I) + C2(J,I) X A2(J,I) = A1(J)*D2(J,I) X 2600 CONTINUE X 2800 CONTINUE X RETURN X END X SUBROUTINE SUB004(M,N,A1,A2,B2,C2,D2) X DIMENSION A2(M,N) X REAL A1(M) X REAL B2(M,N) X REAL C2(M,N) X REAL D2(M,N) X DATA S1/0/ C C DOUBLY DIMENSIONED ARRAY - SIMPLE ASSIGNMENT C X DO 1200 I = 1, N X DO 1000 J = 1, M X A2(J,I) = B2(J,I) + C2(J,I) X 1000 CONTINUE X 1200 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH SCALAR MULTIPLIER C X DO 1600 I = 1, N X DO 1400 J = 1, M X A2(J,I) = B2(J,I) + S1*C2(J,I) X 1400 CONTINUE X 1600 CONTINUE C C DOUBLY DIMENSIONED ARRAY - WITH REUSED VALUE C X DO 2000 I = 1, N X DO 1800 J = 1, M X A2(J,I) = B2(J,I) + C2(J,I) X B2(J,I) = A2(J,I) + D2(J,I) X 1800 CONTINUE X 2000 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH PRIVATE SCALAR TEMPORARY C X DO 2400 I = 1, N X DO 2200 J = 1, M X A06 = B2(J,I) + C2(J,I) X A2(J,I) = A06*D2(J,I) X 2200 CONTINUE X 2400 CONTINUE C C DOUBLY DIMENSIONED ARRAY WITH PRIVATE ARRAY TEMPORARY C X DO 2800 I = 1, N X DO 2600 J = 1, M X A1(J) = B2(J,I) + C2(J,I) X A2(J,I) = A1(J)*D2(J,I) X 2600 CONTINUE X 2800 CONTINUE X RETURN X END X SUBROUTINE SUB005(SUM) X PARAMETER(II=101) X COMMON /DATA/ X 1 R(II,II),Z(II,II),U(II,II),V(II,II),AJ(II,II) X 2,ENERGY(II,II),P(II,II),Q(II,II),TEMP(II,II) X 3,RHO(II,II),DTAU(II,II),MASS(II,II),NBC(II,II) X REAL MASS X COMMON /DATA/ X 1 A(II,II),B(II,II),CBB(II,II),DBB(II,II) X 2,ODTEMP(II,II),SIG(II,II),CC(II,II) X COMMON /DATA/ SUMKV(II),SUMLV(II),DENOMV(II),TEMPRV(II) X DATA TFLR /0.0001/ X LMNP = 2 X KMNP = 2 X LMX = 100 X KMX = 100 X LMN = 2 X KMN = 2 X DTNPH = 3.14 X DO 1200 L = LMNP, LMX X DO 1000 K = KMNP, KMX X CC(K,L) = (0.0001*SQRT(TEMP(K,L))*TEMP(K,L)**2)/AJ(K,L) X SIG(K,L) = MASS(K,L)*SIG(K,L)/DTNPH X ODTEMP(K,L) = TEMP(K,L) X 1000 CONTINUE X 1200 CONTINUE X X X DO 1600 L = LMNP, LMX X DO 1400 K = KMN, KMX X DBB(K,L) = (2.0*CC(K+1,L)*CC(K,L))/(CC(K+1,L)+CC(K,L))*(0.5*(R X 1 (K,L-1)+R(K,L))*((R(K,L)-R(K,L-1))**2+(Z(K,L)-Z(K,L-1))**2) X 2 ) X 1400 CONTINUE X 1600 CONTINUE X X DO 2000 K = KMNP, KMX X DO 1800 L = LMN, LMX X CBB(K,L) = (2.0*CC(K,L)*CC(K,L+1))/(CC(K,L)+CC(K,L+1))*(0.5*(R X 1 (K-1,L)+R(K,L))*((R(K,L)-R(K-1,L))**2+(Z(K,L)-Z(K-1,L))**2) X 2 ) X 1800 CONTINUE X 2000 CONTINUE X X X DO 2400 L = LMNP, LMX X DO 2200 K = KMNP, KMX X DENOMV(K) = SIG(K,L) + CBB(K,L) + CBB(K,L-1)*(1.-A(K,L-1)) X A(K,L) = CBB(K,L)/DENOMV(K) X B(K,L) = (SIG(K,L)*TEMP(K,L)+CBB(K,L-1)*B(K,L-1))/DENOMV(K) X 2200 CONTINUE X 2400 CONTINUE X X ML = LMX + 1 X DO 2800 L = LMNP, LMX X ML = ML - 1 X DO 2600 K = KMNP, KMX X TEMP(K,ML) = A(K,ML)*TEMP(K,ML+1) + B(K,ML) X 2600 CONTINUE X 2800 CONTINUE X X DO 3200 K = KMNP, KMX X DO 3000 L = LMNP, LMX X DENOMV(L) = SIG(K,L) + DBB(K,L) + DBB(K-1,L)*(1.-A(K-1,L)) X A(K,L) = DBB(K,L)/DENOMV(L) X B(K,L) = (SIG(K,L)*TEMP(K,L)+DBB(K-1,L)*B(K-1,L))/DENOMV(L) X 3000 CONTINUE X 3200 CONTINUE X X ML = KMX + 1 X DO 3600 K = KMNP, KMX X ML = ML - 1 X DO 3400 L = LMNP, LMX X TEMP(ML,L) = A(ML,L)*TEMP(ML+1,L) + B(ML,L) X 3400 CONTINUE X 3600 CONTINUE X X YE = -1.0 X DO 4400 L = LMNP, LMX X DO 3800 K = KMNP, KMX X TEMP(K,L) = AMAX1(TEMP(K,L),TFLR) X TEMPRV(K) = ABS((TEMP(K,L)-ODTEMP(K,L))/ODTEMP(K,L)) X 3800 CONTINUE X X DO 4200 K = KMNP, KMX X IF (TEMPRV(K) .LE. YE) GO TO 4000 X YE = TEMPRV(K) X KE = K X LE = L X 4000 CONTINUE X 4200 CONTINUE X X 4400 CONTINUE X KEN = KE X LEN = LE X X SUM = 0.0 X DO 4600 K = KMN, KMX X SUMKV(K) = CBB(K,LMN)*(TEMP(K,LMN)-TEMP(K,LMN+1)) + CBB(K,LMX)*( X 1 TEMP(K,LMX+1)-TEMP(K,LMX)) X 4600 CONTINUE X X DO 4800 L = LMN, LMX X SUMLV(L) = DBB(KMN,L)*(TEMP(KMN,L)-TEMP(KMN+1,L)) + DBB(KMX,L)*( X 1 TEMP(KMX+1,L)-TEMP(KMX,L)) X 4800 CONTINUE X X DO 5000 K = KMN, KMX X SUM = SUM + SUMKV(K) X 5000 CONTINUE X DO 5200 L = LMN, LMX X SUM = SUM + SUMLV(L) X 5200 CONTINUE X X RETURN X END X SUBROUTINE SUB006 X PARAMETER (NY=10,NX=100,NZ=100) X COMMON /DATA/ X * U(NZ),T(NZ), X * UX(NX,NY,NZ),VY(NX,NY,NZ),POTT(NX,NY,NZ),DKZM(NX,NY,NZ), X * DKZH(NX,NY,NZ),WZ(NX,NY,NZ),OBUK(NX,NY),USTR(NX,NY), X * VDEP(NX,NY),TAVR(NX,NY),TSTR(NX,NY),STEPH(NX,NY),Z0(NX,NY), X * ELEV(NX,NY),HMIX(NX,NY),TM(NZ),VM(NZ),UM(NZ) X DO 2000 J = 1, NY X DO 1800 I = 1, NX X ZMH = HMIX(I,J) X ZNOT = Z0(I,J) X IF (ELEV(I,J) .LT. 0) THEN X USTAR = MAX(USTR(I,J),0.001) X ZNOT = 3.905E-5/USTAR + 1.6046E-3*USTAR**2 - 1.747E-4 X Z0(I,J) = ZNOT X ENDIF X IF (J.EQ.1 .OR. J.EQ.NY) THEN X IF (I.EQ.1 .OR. I.EQ.NX) THEN X DO 1000 K = 1, NZ X TX = 0.5*POTT(I,J,K) X WUX = 0.5*UX(I,J,K) X WVX = 0.5*VY(I,J,K) X X TY = 0.5*POTT(I,J,K) X WUY = 0.5*UX(I,J,K) X WVY = 0.5*VY(I,J,K) X XW = (UX(I,J,K)+WUX+WUY)/2. + UM(K) X YW = (VY(I,J,K)+WVX+WVY)/2. + VM(K) X U(K) = SQRT(XW*XW+YW*YW) X T(K) = (POTT(I,J,K)+TX+TY)/2. + TM(K) X 1000 CONTINUE X ELSE X DO 1200 K = 1, NZ X TX = 0.25*(POTT(I-1,J,K)+POTT(I+1,J,K)) X WUX = 0.25*(UX(I-1,J,K)+UX(I+1,J,K)) X WVX = 0.25*(VY(I-1,J,K)+VY(I+1,J,K)) X X TY = 0.5*POTT(I,J,K) X WUY = 0.5*UX(I,J,K) X WVY = 0.5*VY(I,J,K) X XW = (UX(I,J,K)+WUX+WUY)/2. + UM(K) X YW = (VY(I,J,K)+WVX+WVY)/2. + VM(K) X U(K) = SQRT(XW*XW+YW*YW) X T(K) = (POTT(I,J,K)+TX+TY)/2. + TM(K) X 1200 CONTINUE X ENDIF X ELSE X IF (I.EQ.1 .OR. I.EQ.NX) THEN X DO 1400 K = 1, NZ X TX = 0.5*POTT(I,J,K) X WUX = 0.5*UX(I,J,K) X WVX = 0.5*VY(I,J,K) X X TY = 0.25*(POTT(I,J-1,K)+POTT(I,J+1,K)) X WUY = 0.25*(UX(I,J-1,K)+UX(I,J+1,K)) X WVY = 0.25*(VY(I,J-1,K)+VY(I,J+1,K)) X XW = (UX(I,J,K)+WUX+WUY)/2. + UM(K) X YW = (VY(I,J,K)+WVX+WVY)/2. + VM(K) X U(K) = SQRT(XW*XW+YW*YW) X T(K) = (POTT(I,J,K)+TX+TY)/2. + TM(K) X 1400 CONTINUE X ELSE X DO 1600 K = 1, NZ X TX = 0.25*(POTT(I-1,J,K)+POTT(I+1,J,K)) X WUX = 0.25*(UX(I-1,J,K)+UX(I+1,J,K)) X WVX = 0.25*(VY(I-1,J,K)+VY(I+1,J,K)) X X TY = 0.25*(POTT(I,J-1,K)+POTT(I,J+1,K)) X WUY = 0.25*(UX(I,J-1,K)+UX(I,J+1,K)) X WVY = 0.25*(VY(I,J-1,K)+VY(I,J+1,K)) X XW = (UX(I,J,K)+WUX+WUY)/2. + UM(K) X YW = (VY(I,J,K)+WVX+WVY)/2. + VM(K) X U(K) = SQRT(XW*XW+YW*YW) X T(K) = (POTT(I,J,K)+TX+TY)/2. + TM(K) X 1600 CONTINUE X ENDIF X ENDIF X 1800 CONTINUE X 2000 CONTINUE X DO 2600 J = 1, NY X DO 2400 I = 1, NX X DL = OBUK(I,J) X TMES = TAVR(I,J) X USTAR = USTR(I,J) X ZMH = MAX(40.,HMIX(I,J)) X ZMH = MIN(ZMH,2000.) X HMIX(I,J) = ZMH X ZNOT = Z0(I,J) X DEPFAC = USTAR/500. X IF (DL .GE. 0) THEN X VDEP(I,J) = DEPFAC X ELSE X ZC = ZMH/DL X IF (ZC .LT. (-30.)) THEN X VDEP(I,J) = 0.5*DEPFAC*(-ZC)**0.6666667 X ELSE X VDEP(I,J) = DEPFAC*(1.+(-300/DL)**0.6666667) X ENDIF X ENDIF X VDEP(I,J) = MIN(0.004,VDEP(I,J)) X DO 2200 K = 1, NZ X DKZM(I,J,K) = U(K) X DKZH(I,J,K) = T(K) X 2200 CONTINUE X 2400 CONTINUE X 2600 CONTINUE X RETURN X END X SUBROUTINE SUB007 X PARAMETER(II=101) X COMMON /DATA/ X * AJ(II,II), X 1 R(II,II),Z(II,II),U(II,II),V(II,II) X 2,ENERGY(II,II),P(II,II),Q(II,II),TEMP(II,II) X 3,RHO(II,II),DTAU(II,II),MASS(II,II),NBC(II,II) X REAL MASS X X COMMON /DATA/ AU(II),AW(II),AUW(II),AJ1(II),AJ3(II) X 1,VN(II),VNP(II),VOL(II) X DATA P1D6 /0.166666666666667/ X DATA VCUT /1.0E-10/ X KMNP = 2 X LMNP = 2 X LMN = 2 X KMN = 2 X LMX = 100 X KMX = 100 X DTN = .05 X DTNPH = 3.14 X X DO 1200 L = LMN, LMX X DO 1000 K = KMN, KMX X AU(K) = (P(K,L)+Q(K,L))*(Z(K,L-1)-Z(K-1,L)) + (P(K+1,L)+Q(K+1, X 1 L))*(Z(K+1,L)-Z(K,L-1)) + (P(K,L+1)+Q(K,L+1))*(Z(K-1,L)-Z(K X 2 ,L+1)) + (P(K+1,L+1)+Q(K+1,L+1))*(Z(K,L+1)-Z(K+1,L)) X AW(K) = (P(K,L)+Q(K,L))*(R(K,L-1)-R(K-1,L)) + (P(K+1,L)+Q(K+1, X 1 L))*(R(K+1,L)-R(K,L-1)) + (P(K,L+1)+Q(K,L+1))*(R(K-1,L)-R(K X 2 ,L+1)) + (P(K+1,L+1)+Q(K+1,L+1))*(R(K,L+1)-R(K+1,L)) X AUW(K) = RHO(K,L)*AJ(K,L) + RHO(K+1,L)*AJ(K+1,L) + RHO(K,L+1)* X 1 AJ(K,L+1) + RHO(K+1,L+1)*AJ(K+1,L+1) X AUW(K) = 2./AUW(K) X AU(K) = -AU(K)*AUW(K) X AW(K) = AW(K)*AUW(K) X U(K,L) = U(K,L) + DTN*AU(K) X V(K,L) = V(K,L) + DTN*AW(K) X 1000 CONTINUE X 1200 CONTINUE X X DO 1600 L = LMN, LMX X DO 1400 K = KMN, KMX X R(K,L) = R(K,L) + DTNPH*U(K,L) X Z(K,L) = Z(K,L) + DTNPH*V(K,L) X 1400 CONTINUE X 1600 CONTINUE X X DO 2000 L = LMNP, LMX X DO 1800 K = KMNP, KMX X AJ1(K) = R(K,L)*(Z(K-1,L)-Z(K,L-1)) + R(K-1,L)*(Z(K,L-1)-Z(K,L X 1 )) + R(K,L-1)*(Z(K,L)-Z(K-1,L)) X IF (AJ1(K) .EQ. 0) AJ1(K) = P1D6 X AJ3(K) = R(K-1,L)*(Z(K-1,L-1)-Z(K,L-1)) + R(K-1,L-1)*(Z(K,L-1) X 1 -Z(K-1,L)) + R(K,L-1)*(Z(K-1,L)-Z(K-1,L-1)) X IF (AJ3(K) .EQ. 0) AJ3(K) = P1D6 X AJ(K,L) = 0.5*(AJ1(K)+AJ3(K)) X VOL(K) = P1D6*((R(K,L)+R(K-1,L)+R(K,L-1))*AJ1(K)+(R(K-1,L)+R(K X 1 -1,L-1)+R(K,L-1))*AJ3(K)) X VN(K) = 1.0/RHO(K,L) X RHO(K,L) = MASS(K,L)/VOL(K) X VNP(K) = 1.0/RHO(K,L) X DTAU(K,L) = VNP(K) - VN(K) X 1800 CONTINUE X 2000 CONTINUE X X RETURN X END X SUBROUTINE SUB008(NM,N,D,E,Z,IERR) C X DIMENSION D(N),E(N),Z(NM,N) X IERR = 0 C X DO 1000 I = 2, N X E(I-1) = E(I) X 1000 CONTINUE C X F = 0.0 X TST1 = 0.0 X E(N) = 0.0 C X DO 2600 L = 1, N X H = ABS(D(L)) + ABS(E(L)) X TST1 = AMAX1(H,TST1) C X DO 1200 M = L, N X TST2 = TST1 + ABS(E(M)) X IF (TST2 .EQ. TST1) GO TO 1400 X 1200 CONTINUE C X 1400 CONTINUE X IF (M .EQ. L) GO TO 2400 X L1 = L + 1 X L2 = L1 + 1 X G = D(L) X P = (D(L1)-G)/(2.0*E(L)) X R = 3.14 X D(L) = E(L)/(P+SIGN(R,P)) X D(L1) = E(L)*(P+SIGN(R,P)) X DL1 = D(L1) X H = G - D(L) X IF (L2 .GT. N) GO TO 1800 C X DO 1600 I = L2, N X D(I) = D(I) - H X 1600 CONTINUE C X 1800 CONTINUE X F = F + H X P = D(M) X C = 1.0 X C2 = C X EL1 = E(L1) X S = 0.0 X MML = M - L X DO 2200 II = 1, MML X C3 = C2 X C2 = C X S2 = S X I = M - II X G = C*E(I) X H = C*P X R = 3.14 X E(I+1) = S*R X S = E(I)/R X C = P/R X P = C*D(I) - S*G X D(I+1) = H + S*(C*G+S*D(I)) X DO 2000 K = 1, N X H = Z(K,I+1) X Z(K,I+1) = S*Z(K,I) + C*H X Z(K,I) = C*Z(K,I) - S*H X 2000 CONTINUE C X 2200 CONTINUE C X P = -S*S2*C3*EL1*E(L)/DL1 X E(L) = S*P X D(L) = C*P X TST2 = TST1 + ABS(E(L)) X 2400 CONTINUE X D(L) = D(L) + F X 2600 CONTINUE X X DO 3600 II = 2, N X I = II - 1 X K = I X P = D(I) C X DO 3000 J = II, N X IF (D(J) .GE. P) GO TO 2800 X K = J X P = D(J) X 2800 CONTINUE X 3000 CONTINUE C X IF (K .EQ. I) GO TO 3400 X D(K) = D(I) X D(I) = P C X DO 3200 J = 1, N X P = Z(J,I) X Z(J,I) = Z(J,K) X Z(J,K) = P X 3200 CONTINUE C X 3400 CONTINUE X 3600 CONTINUE C X RETURN X END X SUBROUTINE SUB009(S,T,U,FNDS) X COMMON /DATA/ FN(27000),SQ(30),TQ(30),UQ(30) X 1 , DSQ(30),DTQ(30),DUQ(30) X DIMENSION FNDS(27000,3) X X A = 0.5 X B = S*S X SQ(1) = A*(B-S) X SQ(2) = A + A - B X SQ(3) = A*(B+S) X B = T*T X TQ(1) = A*(B-T) X TQ(2) = A + A - B X TQ(3) = A*(B+T) X B = U*U X UQ(1) = A*(B-U) X UQ(2) = A + A - B X UQ(3) = A*(B+U) X X L = 0 X DO 1400 IU = 1, 30 X UQI = UQ(IU) X DO 1200 IT = 1, 30 X TQI = TQ(IT) X DO 1000 IS = 1, 30 X L = L + 1 X FN(L) = UQI*TQI*SQ(IS) X 1000 CONTINUE X 1200 CONTINUE X 1400 CONTINUE X X A = 0.5 X DSQ(1) = S - A X DSQ(2) = (-S) - S X DSQ(3) = S + A X DTQ(1) = T - A X DTQ(2) = (-T) - T X DTQ(3) = T + A X DUQ(1) = U - A X DUQ(2) = (-U) - U X DUQ(3) = U + A X X L = 0 X DO 2000 IU = 1, 30 X UQI = UQ(IU) X DUQI = DUQ(IU) X DO 1800 IT = 1, 30 X TQI = TQ(IT) X DTQI = DTQ(IT) X DO 1600 IS = 1, 30 X L = L + 1 X FNDS(L,1) = UQI*TQI*DSQ(IS) X FNDS(L,2) = UQI*DTQI*SQ(IS) X FNDS(L,3) = DUQI*TQI*SQ(IS) X 1600 CONTINUE X 1800 CONTINUE X 2000 CONTINUE X END X SUBROUTINE SUB010(NSRC,IPAR,ISLECT) X PARAMETER (ISIZE1=12, ISIZE2=24) X PARAMETER (MAXVOL=ISIZE1**3*ISIZE2, NC=3, NC8=8*NC) X PARAMETER (NN1=ISIZE1/2, NN48=NC8*MAXVOL/2, NC4=4*NC) X PARAMETER (NSC=ISIZE1**3) X COMMON /DATA/ XX(NC8,NN1,ISIZE1,ISIZE1,ISIZE2) X COMMON /DATA/ PP(NC8,NN1,ISIZE1,ISIZE1,ISIZE2) X COMMON /DATA/ RQ(NN48) X COMMON /DATA/ PH(2,NC,2,2,NSC) X COMMON /DATA/ PN(2,NC,2,2,NSC) X X DIMENSION PNS(NC8,ISIZE1,ISIZE1,ISIZE1) X DIMENSION ISLECT(4),RX(NN48) X DIMENSION PQ(2,NC,2,2,NSC,ISIZE2) X EQUIVALENCE (RX,XX), (PP,PQ) X EQUIVALENCE(PN,PNS) X B0 = 3.14 X DO 6600 IT = 1, ISIZE2 X IF (NSRC .EQ. 1) THEN X DO 1800 IMRE = 1, 2 X DO 1600 MUS = 1, 2 X DO 1400 MUG = 1, 2 X DO 1200 IC = 1, NC X DO 1000 ISC = 1, NSC X PH(IMRE,IC,MUS,MUG,ISC) = PQ(IMRE,IC,MUS,MUG,ISC, X 1 IT) X 1000 CONTINUE X 1200 CONTINUE X 1400 CONTINUE X 1600 CONTINUE X 1800 CONTINUE X ELSE IF (NSRC .EQ. 5) THEN X DO 2400 IMRE = 1, 2 X DO 2200 IC = 1, 3 X DO 2000 ISC = 1, NSC X PH(IMRE,IC,1,1,ISC) = PQ(IMRE,IC,1,2,ISC,IT) X PH(IMRE,IC,2,1,ISC) = PQ(IMRE,IC,2,2,ISC,IT) X PH(IMRE,IC,1,2,ISC) = PQ(IMRE,IC,1,1,ISC,IT) X PH(IMRE,IC,2,2,ISC) = PQ(IMRE,IC,2,1,ISC,IT) X 2000 CONTINUE X 2200 CONTINUE X 2400 CONTINUE X ENDIF X IF (IT.NE.1 .AND. ISELC.EQ.0) GO TO 5400 X DO 3400 IMRE = 1, 2 X DO 3200 IC = 1, NC X DO 3000 MUS = 1, 2 X DO 2800 MUG = 1, 2 X DO 2600 ISC = 1, NSC X PN(IMRE,IC,MUS,MUG,ISC) = 0. X 2600 CONTINUE X 2800 CONTINUE X 3000 CONTINUE X 3200 CONTINUE X 3400 CONTINUE X ISELC = 0 X DO 5200 IS = 1, 4 X IF (ISLECT(IS) .EQ. 1) THEN X DO 3800 IC = 1, NC X DO 3600 ISC = 1, NSC X PN(1,IC,1,1,ISC) = PN(1,IC,1,1,ISC) + PH(2,IC,2,2,ISC) X PN(2,IC,1,1,ISC) = PN(2,IC,1,1,ISC) - PH(1,IC,2,2,ISC) X PN(1,IC,2,1,ISC) = PN(1,IC,2,1,ISC) + PH(2,IC,1,2,ISC) X PN(2,IC,2,1,ISC) = PN(2,IC,2,1,ISC) - PH(1,IC,1,2,ISC) X PN(1,IC,1,2,ISC) = PN(1,IC,1,2,ISC) - PH(2,IC,2,1,ISC) X PN(2,IC,1,2,ISC) = PN(2,IC,1,2,ISC) + PH(1,IC,2,1,ISC) X PN(1,IC,2,2,ISC) = PN(1,IC,2,2,ISC) - PH(2,IC,1,1,ISC) X PN(2,IC,2,2,ISC) = PN(2,IC,2,2,ISC) + PH(1,IC,1,1,ISC) X 3600 CONTINUE X 3800 CONTINUE X ISELC = ISELC + 1 X ELSE IF (ISLECT(IS) .EQ. 2) THEN X DO 4200 IC = 1, NC X DO 4000 ISC = 1, NSC X PN(1,IC,1,1,ISC) = PN(1,IC,1,1,ISC) - PH(1,IC,2,2,ISC) X PN(2,IC,1,1,ISC) = PN(2,IC,1,1,ISC) - PH(2,IC,2,2,ISC) X PN(1,IC,2,1,ISC) = PN(1,IC,2,1,ISC) + PH(1,IC,1,2,ISC) X PN(2,IC,2,1,ISC) = PN(2,IC,2,1,ISC) + PH(2,IC,1,2,ISC) X PN(1,IC,1,2,ISC) = PN(1,IC,1,2,ISC) + PH(1,IC,2,1,ISC) X PN(2,IC,1,2,ISC) = PN(2,IC,1,2,ISC) + PH(2,IC,2,1,ISC) X PN(1,IC,2,2,ISC) = PN(1,IC,2,2,ISC) - PH(1,IC,1,1,ISC) X PN(2,IC,2,2,ISC) = PN(2,IC,2,2,ISC) - PH(2,IC,1,1,ISC) X 4000 CONTINUE X 4200 CONTINUE X ISELC = ISELC + 1 X ELSE IF (ISLECT(IS) .EQ. 3) THEN X DO 4600 IC = 1, NC X DO 4400 ISC = 1, NSC X PN(1,IC,1,1,ISC) = PN(1,IC,1,1,ISC) + PH(2,IC,1,2,ISC) X PN(2,IC,1,1,ISC) = PN(2,IC,1,1,ISC) - PH(1,IC,1,2,ISC) X PN(1,IC,2,1,ISC) = PN(1,IC,2,1,ISC) - PH(2,IC,2,2,ISC) X PN(2,IC,2,1,ISC) = PN(2,IC,2,1,ISC) + PH(1,IC,2,2,ISC) X PN(1,IC,1,2,ISC) = PN(1,IC,1,2,ISC) - PH(2,IC,1,1,ISC) X PN(2,IC,1,2,ISC) = PN(2,IC,1,2,ISC) + PH(1,IC,1,1,ISC) X PN(1,IC,2,2,ISC) = PN(1,IC,2,2,ISC) + PH(2,IC,2,1,ISC) X PN(2,IC,2,2,ISC) = PN(2,IC,2,2,ISC) - PH(1,IC,2,1,ISC) X 4400 CONTINUE X 4600 CONTINUE X ISELC = ISELC + 1 X ELSE IF (ISLECT(IS) .EQ. 4) THEN X DO 5000 IC = 1, NC X DO 4800 ISC = 1, NSC X PN(1,IC,1,1,ISC) = PN(1,IC,1,1,ISC) - PH(1,IC,1,1,ISC) X PN(2,IC,1,1,ISC) = PN(2,IC,1,1,ISC) - PH(2,IC,1,1,ISC) X PN(1,IC,2,1,ISC) = PN(1,IC,2,1,ISC) - PH(1,IC,2,1,ISC) X PN(2,IC,2,1,ISC) = PN(2,IC,2,1,ISC) - PH(2,IC,2,1,ISC) X PN(1,IC,1,2,ISC) = PN(1,IC,1,2,ISC) + PH(1,IC,1,2,ISC) X PN(2,IC,1,2,ISC) = PN(2,IC,1,2,ISC) + PH(2,IC,1,2,ISC) X PN(1,IC,2,2,ISC) = PN(1,IC,2,2,ISC) + PH(1,IC,2,2,ISC) X PN(2,IC,2,2,ISC) = PN(2,IC,2,2,ISC) + PH(2,IC,2,2,ISC) X 4800 CONTINUE X 5000 CONTINUE X ISELC = ISELC + 1 X ENDIF X 5200 CONTINUE X 5400 CONTINUE X IF (ISELC .EQ. 0) THEN X DO 5600 IALL = 1, 8*NC*NSC X PN(IALL,1,1,1,1) = PH(IALL,1,1,1,1) X 5600 CONTINUE X ENDIF X DO 6400 IZ = 1, ISIZE1 X DO 6200 IY = 1, ISIZE1 X DO 6000 IX = 1, NN1 X J0 = MOD(IY+IZ+IT+IPAR,2) X J1 = 1 - J0 X DO 5800 I = 1, NC8 X XX(I,IX,IY,IZ,IT) = PNS(I,2*IX-J0,IY,IZ) X PP(I,IX,IY,IZ,IT) = PNS(I,2*IX-J1,IY,IZ) X 5800 CONTINUE X 6000 CONTINUE X 6200 CONTINUE X 6400 CONTINUE X 6600 CONTINUE X DO 6800 I = 1, NN48 X RX(I) = B0*RX(I) - RQ(I) X 6800 CONTINUE X RETURN X END X SUBROUTINE SUB011 X PARAMETER (N=1000) X COMMON /DATA/ B(N), A(N,N) X DO 1200 I = 2, N X SUM = 0. X DO 1000 L = 1, I - 1 X SUM = SUM + A(L,I)*B(L) X 1000 CONTINUE X B(I) = B(I) - SUM X 1200 CONTINUE X X NRMAND = (N-1)*1 X IF (NRMAND .EQ. 0) THEN X X DO 1600 J = 1, N - 1, 2 X B(J+1) = B(J+1) - A(J,J+1)*B(J) X DO 1400 I = J + 2, N X B(I) = B(I) - A(J,I)*B(J) - A(J+1,I)*B(J+1) X 1400 CONTINUE X 1600 CONTINUE X ELSE X X DO 1800 I = 2, N X B(I) = B(I) - A(1,I)*B(1) X 1800 CONTINUE X X DO 2200 J = 2, N - 1, 2 X B(J+1) = B(J+1) - A(J,J+1)*B(J) X DO 2000 I = J + 2, N X B(I) = B(I) - A(J,I)*B(J) - A(J+1,I)*B(J+1) X 2000 CONTINUE X 2200 CONTINUE X X ENDIF X X K = 1 X DO 2400 I = 1, N X B(I) = B(I)*A(I,I) X B(I) = B(I)*A(K,1) X K = K + N + 1 X 2400 CONTINUE X DO 2800 I = N - 1, 1, -1 X SUM = 0. X DO 2600 L = I + 1, N X SUM = SUM + A(I,L)*B(L) X 2600 CONTINUE X B(I) = B(I) - SUM X 2800 CONTINUE C X IF (NRMAND .EQ. 0) THEN X X DO 3200 J = N, 2, -2 X B(J-1) = B(J-1) - A(J-1,J)*B(J) X DO 3000 I = J - 2, 1, -1 X B(I) = B(I) - A(I,J)*B(J) - A(I,J-1)*B(J-1) X 3000 CONTINUE X 3200 CONTINUE X X ELSE X X DO 3400 I = N - 1, 1, -1 X B(I) = B(I) - A(I,N)*B(N) X 3400 CONTINUE X X DO 3800 J = N - 1, 2, -2 X B(J-1) = B(J-1) - A(J-1,J)*B(J) X DO 3600 I = J - 2, 1, -1 X B(I) = B(I) - A(I,J)*B(J) - A(I,J-1)*B(J-1) X 3600 CONTINUE X 3800 CONTINUE X X ENDIF X X X RETURN X END X SUBROUTINE SUB012(B) X PARAMETER (IDA=100,M=100,N=25,NRHS=100,IDB=100,NMAT=100) X COMMON/DATA/ A(0:IDA, -M:0, 0:N),EPSS(0:256) X DIMENSION B(0:NRHS, 0:IDB, 0:N) X DATA EPS/1E-13/ X X DO 2600 J = 0, N X I0 = MAX((-M),(-J)) X X DO 1600 I = I0, -1 X DO 1200 JJ = I0 - I, -1 X DO 1000 L = 0, NMAT X A(L,I,J) = A(L,I,J) - A(L,JJ,I+J)*A(L,I+JJ,J) X 1000 CONTINUE X 1200 CONTINUE X DO 1400 L = 0, NMAT X A(L,I,J) = A(L,I,J)*A(L,0,I+J) X 1400 CONTINUE X 1600 CONTINUE X X DO 1800 L = 0, NMAT X EPSS(L) = EPS*A(L,0,J) X 1800 CONTINUE X DO 2200 JJ = I0, -1 X DO 2000 L = 0, NMAT X A(L,0,J) = A(L,0,J) - A(L,JJ,J)**2 X 2000 CONTINUE X 2200 CONTINUE X DO 2400 L = 0, NMAT X A(L,0,J) = 1./SQRT(ABS(EPSS(L)+A(L,0,J))) X 2400 CONTINUE X 2600 CONTINUE X X DO 4400 I = 0, NRHS X DO 3400 K = 0, N X DO 2800 L = 0, NMAT X B(I,L,K) = B(I,L,K)*A(L,0,K) X 2800 CONTINUE X DO 3200 JJ = 1, MIN(M,N-K) X DO 3000 L = 0, NMAT X B(I,L,K+JJ) = B(I,L,K+JJ) - A(L,(-JJ),K+JJ)*B(I,L,K) X 3000 CONTINUE X 3200 CONTINUE X 3400 CONTINUE X X DO 4200 K = N, 0, -1 X DO 3600 L = 0, NMAT X B(I,L,K) = B(I,L,K)*A(L,0,K) X 3600 CONTINUE X DO 4000 JJ = 1, MIN(M,K) X DO 3800 L = 0, NMAT X B(I,L,K-JJ) = B(I,L,K-JJ) - A(L,(-JJ),K)*B(I,L,K) X 3800 CONTINUE X 4000 CONTINUE X 4200 CONTINUE X 4400 CONTINUE X X RETURN X END X SUBROUTINE SUB013 X PARAMETER (IMAX=1000,KMAX=100) X COMMON /DATA/ GESHEM( IMAX ) X COMMON /DATA/ PS( IMAX ),DEL( IMAX ),SL( IMAX ), X 1QN( IMAX , KMAX ),QN1( IMAX , KMAX ),DQ( IMAX , KMAX ), X 2TN1( IMAX , KMAX ) X COMMON /DATA/ ACUM( IMAX ),PRESS( KMAX ),TIN( KMAX ),QIN( KMAX ), X 1TMST( KMAX ),QMST( KMAX ),DTKUO( KMAX ),DQKUO( KMAX ) X 2,ESAT( KMAX ) X X DT = 3.14 X ACUMT = 4.15 X NKUO = 0 X RELKUO = .2175 X MSTA = 0 X RDT = 1./DT X CPOVL = 1005./2.5E+6 X RELEPS = RELKUO*0.622 X DO 1200 K = 1, KMAX X DO 1000 I = 1, IMAX X DQ(I,K) = QN1(I,K) - QN(I,K) X 1000 CONTINUE X 1200 CONTINUE X KACUM = 4 X DO 1400 I = 1, IMAX X ACUM(I) = ACUMT*0.1 X 1400 CONTINUE X DO 1800 K = 1, KMAX X DO 1600 I = 1, IMAX X ACUM(I) = ACUM(I) + RDT*DQ(I,K)*DEL(K) X 1600 CONTINUE X 1800 CONTINUE X DO 4200 I = 1, IMAX X DO 2000 K = 1, KMAX X PRESS(K) = SL(K)*PS(I) X 2000 CONTINUE X DO 2200 K = 1, KMAX X TIN(K) = TN1(I,K) X QIN(K) = QN1(I,K) X IF (QN1(I,K) .LT. .0) QIN(K) = .0 X 2200 CONTINUE X QSATK = .65*.622*ESAT(1)/(PRESS(1)-0.378*ESAT(1)) X QDEF = QIN(1) - QSATK X TIN(1) = TIN(1) + .25*(TIN(2)-TIN(1)) X HNEW = PS(I)/100. X DO 2400 K = 2, KMAX X X = TMST(K) - TIN(K) X 2400 CONTINUE X NKUO = NKUO + 1 X WATER = 0. X WATER = WATER + DQ(I,1)*DEL(1) X WATER = WATER + DQ(I,2)*DEL(2) X Q1 = 0. X Q2 = 0. X DO 3000 K = 2, KMAX X X = QMST(K) - QIN(K) X IF (X .LE. 0) GO TO 2600 X Q1 = Q1 + X*DEL(K) X 2600 CONTINUE X X = TMST(K) - TIN(K) X IF (X .LE. 0) GO TO 2800 X Q2 = Q2 + X*DEL(K) X 2800 CONTINUE X 3000 CONTINUE X Q2 = Q2*CPOVL X QEFF = .92 X IF (Q1.GT.0 .OR. Q2.GT.0) QEFF = WATER/(Q1+Q2) X QEFF = AMIN1(1.00,QEFF) X 3200 CONTINUE X DTKUO(K) = X X X = QEFF*(QMST(K)-QIN(K)) X IF (X .GT. 0) GO TO 3400 X X = 0. X 3400 CONTINUE X DQKUO(K) = X X 3600 CONTINUE X DTKUO(1) = 0. X DQKUO(1) = 0. X CEVAP = 0. X DO 3800 K = 2, KMAX X TIN(K) = TIN(K) + DTKUO(K) X QIN(K) = DQKUO(K) + QIN(K) X CEVAP = CEVAP + DQKUO(K)*DEL(K) X 3800 CONTINUE X FALL = WATER - CEVAP X GESHEM(I) = GESHEM(I) + PS(I)*FALL/9.8 X QIN(1) = QIN(1) - DQ(I,1) X IF (DQ(I,2) .GT. 0.) QIN(2) = QIN(2) - DQ(I,2) X TIN(1) = TN1(I,1) X DO 4000 K = 1, KMAX X QN1(I,K) = QIN(K) X TN1(I,K) = TIN(K) X 4000 CONTINUE X 4200 CONTINUE X RETURN X END X SUBROUTINE SUB014 X PARAMETER (NI=48,NJ=72,NK=32,NJK=NJ*NK,DTH=1./NJ,DPH=1./NK X $ ,DS=1./NI,DTDP=DTH*DPH,DVOL=DS*DTDP,TSI=NI,TDS=.5*TSI X $ ,MLMNV=50,MLMNB=160,NMIN=-4,NMAX=4,MM=18,QN=1.00,QB=0.0 X $ ,NSTA=4,DSTA=2.*DTDP/NSTA,MLMNS=32,MMS=15,INSOL=0 X $ ,IVAC=0,PVAC=1.00,NVI=NI+IVAC X $ ,NSMIN=-9,NSMAX=13,NJKS=NJK*NSTA,NPROCS=1,QONAX=9.99 X $ ,NSKIP=1,FSCALE=1.7,FXIN=5.0,PSHIFT=1.0,PARFAC=0.5 ) X X PARAMETER (MD=MLMNS,MDY=MLMNS,ND=NVI,ND1=ND+1,NA=(MD+MDY)*ND+MD) X X PARAMETER(LMNB=160) X COMMON /DATA/BS(NJK,0:NI),TCOS(NJK,MLMNB),TSIN(NJK,MLMNB) X $ ,VJAC(NJK,0:NI) X COMMON /DATA/ FPHV (MLMNB,NI) X $ ,FBJAC(MLMNB,0:NI), FBS(MLMNB,0:NI) X COMMON /DATA/ BJAC(NJK,0:NI),BJACS(NJK,0:NI) X $ , GTTL(NJK,0:NI), GTPL(NJK,0:NI) X $ , GPPL(NJK,0:NI) X COMMON /DATA/ BT(NJK,0:NI), BP(NJK,0:NI), BSQ(NJK,0:NI) X COMMON /DATA/ AM( NI), PTH( NI), PVP( NI),PVPI( NI) X $ , FPP( NI), FTP( NI), PPI( NI), PP( NI) X $ , FPPP( NI), FTPP( NI),AIOTA( NI), VVP( NI) X $ , WMAGV( NI), WMAG( NI) X $ , CIV( NI), CJV( NI), CIVP( NI),CJVP( NI) X $ , CI ( NI), CJ ( NI), CIP ( NI),CJP ( NI) X $ , CIPI( NI), CJPI( NI) X $ , VP ( NI), VPP( NI),EQUIV( NI),EQUI( NI) X LVMTPR = 1 X LCURR = 1 X DO 1400 I = 1, NI X CI(I) = 0. X CJ(I) = 0. X VP(I) = 0. X WMAG(I) = 0. X DO 1000 JK = 1, NJK X BP(JK,I) = GPPL(JK,I)*FTP(I) + GTPL(JK,I)*FPP(I) X BT(JK,I) = GTPL(JK,I)*FTP(I) + GTTL(JK,I)*FPP(I) X IF (LCURR .EQ. 1) THEN X CI(I) = CI(I) - DTDP*BP(JK,I) X CJ(I) = CJ(I) + DTDP*BT(JK,I) X ENDIF X IF (LCURR .EQ. 2) THEN X CI(I) = CIV(I) X CJ(I) = CJV(I) X ENDIF X FPP(I) = FPP(I) X VP(I) = VP(I) - DTDP*BJAC(JK,I) X 1000 CONTINUE X DO 1200 JK = 1, NJK X WMAG(I) = WMAG(I) + 0.5*BSQ(JK,I)*BJAC(JK,I)*DVOL X 1200 CONTINUE X 1400 CONTINUE X DO 1600 I = 1, NI - 1 X VPI = 0.5*(VP(I+1)+VP(I)) X FTPI = 0.5*(FTP(I+1)+FTP(I)) X FPPI = 0.5*(FPP(I+1)+FPP(I)) X CIPI(I) = TSI*(CI(I+1)-CI(I)) X CJPI(I) = TSI*(CJ(I+1)-CJ(I)) X PVPI(I) = TSI*(PTH(I+1)-PTH(I)) X EQUIN = ABS(CJPI(I)*FPPI) + ABS(CIPI(I)*FTPI) + ABS(VPI*PVPI(I)) X EQUI(I) = ((CJPI(I)*FPPI-CIPI(I)*FTPI)-VPI*PVPI(I))/EQUIN X PPI(I) = (CJPI(I)*FPPI-CIPI(I)*FTPI)/VPI X 1600 CONTINUE X DO 2000 I = 2, NI - 1 X FPPP(I) = TDS*(FPP(I+1)-FPP(I-1)) X FTPP(I) = TDS*(FTP(I+1)-FTP(I-1)) X PP(I) = 0.5*(PPI(I)+PPI(I-1)) X CIP(I) = 0.5*(CIPI(I)+CIPI(I-1)) X CJP(I) = 0.5*(CJPI(I)+CJPI(I-1)) X DO 1800 JK = 1, NJK X BJACS(JK,I) = TDS*(BJAC(JK,I+1)-BJAC(JK,I-1)) X 1800 CONTINUE X 2000 CONTINUE X FTPP(NI) = 2.0*FTPP(NI-1) - FTPP(NI-2) X FPPP(NI) = 2.0*FPPP(NI-1) - FPPP(NI-2) X PP(NI) = 1.5*PPI(NI-1) - PPI(NI-2)*0.5 X CIP(NI) = 1.5*CIPI(NI-1) - CIPI(NI-2)*0.5 X CJP(NI) = 1.5*CJPI(NI-1) - CJPI(NI-2)*0.5 X FTPP(1) = 2.0*FTPP(2) - FTPP(3) X FPPP(1) = 2.0*FPPP(2) - FPPP(3) X PP(1) = 1.5*PPI(1) - PPI(2)*0.5 X CIP(1) = 1.5*CIPI(1) - CIPI(2)*0.5 X CJP(1) = 1.5*CJPI(1) - CJPI(2)*0.5 X DO 2200 JK = 1, NJK X BJACS(JK,1) = 2.0*BJACS(JK,2) - BJACS(JK,3) X BJACS(JK,NI) = 2.0*BJACS(JK,NI-1) - BJACS(JK,NI-2) X 2200 CONTINUE X DO 3600 I = 1, NI X DO 2400 JK = 1, NJK X VJAC(JK,I) = (FPP(I)*CJ(I)-FTP(I)*CI(I))/BSQ(JK,I) X 2400 CONTINUE X DO 2800 L = 1, LMNB X FPHV(L,I) = 0. X FBS(L,I) = 0. X DO 2600 JK = 1, NJK X FPHV(L,I) = FPHV(L,I) + VJAC(JK,I)*TCOS(JK,L) X 2600 CONTINUE X FPHV(L,I) = 2.*DTDP*FPHV(L,I) X DIFJAC = 2.*(FBJAC(L,I)-FPHV(L,I))/(VP(I)+VVP(I)) X IF (ABS(DIFJAC).GE.1.E-06 .AND. LVMTPR.EQ.1) THEN X ENDIF X 2800 CONTINUE X DO 3000 JK = 1, NJK X BS(JK,I) = 0. X 3000 CONTINUE X DO 3400 L = 1, LMNB X DO 3200 JK = 1, NJK X BS(JK,I) = BS(JK,I) + FBS(L,I)*TSIN(JK,L) X 3200 CONTINUE X 3400 CONTINUE X 3600 CONTINUE X DO 3800 JK = 1, NJK X BS(JK,1) = 2.*BS(JK,2) - BS(JK,3) X 3800 CONTINUE X RETURN X END X SUBROUTINE SUB015(FX,JSAVE,KC) X PARAMETER (NOMOL=172,NATOMS=3,M=10000) X PARAMETER (NMOL=172,NATOM=3,NMOL1=171,NATMO=1) X PARAMETER (MAX = NOMOL*NATOM,MXCPUS=8,NCPUS=8) X COMMON /DATA/ TEMP,RHO,TSTEP,BOXL,BOXH,CUTOFF,CUT2 X COMMON /DATA/ OMAS,HMAS,WTMOL,ROH,ANGLE,FHM,FOM,ROHI,ROHI2 X COMMON /DATA/ QQ,A1,B1,A2,B2,A3,B3,A4,B4,AB1,AB2,AB3,AB4,C1,C2, X * QQ2,QQ4,REF1,REF2,REF4 X COMMON /DATA/ UNITT,UNITL,UNITM,BOLTZ,AVGNO,PCC(10) X COMMON /DATA/ X(M),Y(M),Z(M),FY(M),FZ(M),XM(M),YM(M),ZM(M) X COMMON /DATA/ XL(NOMOL,14),YL(NOMOL,14),ZL(NOMOL,14),RL(NOMOL,14), X . RS(NOMOL,14),FF(NOMOL,14),GG(NOMOL,14),FTEMP(NOMOL,4) X COMMON /DATA/ VIRCPU(MXCPUS),FXT(MAX,MXCPUS), X . FYT(MAX,MXCPUS),FZT(MAX,MXCPUS) X DIMENSION FX(M),JSAVE(NOMOL),KC(NOMOL) X FUNEXP(AB,B,XX) = AB*EXP(-B*XX)*XX X IDCPU = 0 X VIR = 0 X DO 1200 ID = 1, NCPUS X VIRCPU(ID) = 0. X DO 1000 I = 1, MAX X FXT(I,ID) = 0. X FYT(I,ID) = 0. X FZT(I,ID) = 0. X 1000 CONTINUE X 1200 CONTINUE X X X IDCPU = IDCPU + 1 X IW10 = 1 X IWO0 = 2 X IW20 = 3 X X DO 6400 I = 1, NMOL1 X X X IW1 = IW10 + (I-1)*NATOMS X IWO = IWO0 + (I-1)*NATOMS X IW2 = IW20 + (I-1)*NATOMS X X DO 1400 J = I + 1, NMOL X X JW1 = IW10 + (J-1)*NATOMS X X XL(J,1) = XM(I) - XM(J) X XL(J,4) = X(IW1) - XM(J) X XL(J,5) = X(IW1+2) - XM(J) X XL(J,2) = XM(I) - X(JW1) X XL(J,6) = X(IW1) - X(JW1) X XL(J,8) = X(IW1+2) - X(JW1) X XL(J,11) = X(IW1+1) - X(JW1) X XL(J,10) = X(IW1+1) - X(JW1+1) X XL(J,13) = X(IW1) - X(JW1+1) X XL(J,14) = X(IW1+2) - X(JW1+1) X XL(J,3) = XM(I) - X(JW1+2) X XL(J,7) = X(IW1) - X(JW1+2) X XL(J,9) = X(IW1+2) - X(JW1+2) X XL(J,12) = X(IW1+1) - X(JW1+2) X X YL(J,1) = YM(I) - YM(J) X YL(J,4) = Y(IW1) - YM(J) X YL(J,5) = Y(IW1+2) - YM(J) X YL(J,2) = YM(I) - Y(JW1) X YL(J,6) = Y(IW1) - Y(JW1) X YL(J,8) = Y(IW1+2) - Y(JW1) X YL(J,11) = Y(IW1+1) - Y(JW1) X YL(J,10) = Y(IW1+1) - Y(JW1+1) X YL(J,13) = Y(IW1) - Y(JW1+1) X YL(J,14) = Y(IW1+2) - Y(JW1+1) X YL(J,3) = YM(I) - Y(JW1+2) X YL(J,7) = Y(IW1) - Y(JW1+2) X YL(J,9) = Y(IW1+2) - Y(JW1+2) X YL(J,12) = Y(IW1+1) - Y(JW1+2) X X ZL(J,1) = ZM(I) - ZM(J) X ZL(J,4) = Z(IW1) - ZM(J) X ZL(J,5) = Z(IW1+2) - ZM(J) X ZL(J,2) = ZM(I) - Z(JW1) X ZL(J,8) = Z(IW1+2) - Z(JW1) X ZL(J,6) = Z(IW1) - Z(JW1) X ZL(J,11) = Z(IW1+1) - Z(JW1) X ZL(J,10) = Z(IW1+1) - Z(JW1+1) X ZL(J,13) = Z(IW1) - Z(JW1+1) X ZL(J,14) = Z(IW1+2) - Z(JW1+1) X ZL(J,3) = ZM(I) - Z(JW1+2) X ZL(J,7) = Z(IW1) - Z(JW1+2) X ZL(J,9) = Z(IW1+2) - Z(JW1+2) X ZL(J,12) = Z(IW1+1) - Z(JW1+2) X 1400 CONTINUE X X DO 1800 K = 1, 14 X DO 1600 J = I + 1, NMOL X IF(ABS(XL(J,K)).GT.BOXH)XL(J,K)=XL(J,K)-SIGN(BOXL,XL(J,K)) X IF(ABS(YL(J,K)).GT.BOXH)YL(J,K)=YL(J,K)-SIGN(BOXL,YL(J,K)) X IF(ABS(ZL(J,K)).GT.BOXH)ZL(J,K)=ZL(J,K)-SIGN(BOXL,ZL(J,K)) X 1600 CONTINUE X 1800 CONTINUE X X DO 2000 J = I + 1, NMOL X KC(J) = 0 X 2000 CONTINUE X X DO 2400 K = 1, 9 X DO 2200 J = I + 1, NMOL X RS(J,K)=XL(J,K)*XL(J,K)+YL(J,K)*YL(J,K)+ZL(J,K)*ZL(J,K) X IF (RS(J,K) .GT. CUT2) KC(J) = KC(J) + 1 X 2200 CONTINUE X 2400 CONTINUE X X JT = 0 X DO 2600 J = I + 1, NMOL X IF (KC(J) .LT. 9) THEN X JT = JT + 1 X JSAVE(JT) = J X ENDIF X 2600 CONTINUE X JNUM = JT X X DO 3000 K = 1, 14 X DO 2800 J = I + 1, NMOL X FF(J,K) = 0. X 2800 CONTINUE X 3000 CONTINUE X X DO 3200 J = I + 1, NMOL X IF(RS(J,1).LT.CUT2)FF(J,1)=QQ4*(RS(J,1)*SQRT(RS(J,1)))+REF4 X IF (RS(J,2) .LT. CUT2) FF(J,2) = (-QQ2*(RS(J,2)*SQRT(RS(J,2))) X 1 ) - REF2 X IF (RS(J,3) .LT. CUT2) FF(J,3) = (-QQ2*(RS(J,3)*SQRT(RS(J,3))) X 1 ) - REF2 X IF (RS(J,4) .LT. CUT2) FF(J,4) = (-QQ2*(RS(J,4)*SQRT(RS(J,4))) X 1 ) - REF2 X IF (RS(J,5) .LT. CUT2) FF(J,5) = (-QQ2*(RS(J,5)*SQRT(RS(J,5))) X 1 ) - REF2 X IF (RS(J,6) .LE. CUT2) THEN X RL(J,6) = SQRT(RS(J,6)) X FF(J,6) = QQ*(RS(J,6)*RL(J,6)) + REF1 X ENDIF X IF (RS(J,7) .LE. CUT2) THEN X RL(J,7) = SQRT(RS(J,7)) X FF(J,7) = QQ*(RS(J,7)*RL(J,7)) + REF1 X ENDIF X IF (RS(J,8) .LE. CUT2) THEN X RL(J,8) = SQRT(RS(J,8)) X FF(J,8) = QQ*(RS(J,8)*RL(J,8)) + REF1 X ENDIF X IF (RS(J,9) .LE. CUT2) THEN X RL(J,9) = SQRT(RS(J,9)) X FF(J,9) = QQ*(RS(J,9)*RL(J,9)) + REF1 X ENDIF X 3200 CONTINUE X X DO 3400 J = I + 1, NMOL X IF (KC(J) .LT. 9) VIRCPU(IDCPU) = VIRCPU(IDCPU) + FF(J,1)*RS(J X 1 ,1) + FF(J,2)*RS(J,2) + FF(J,3)*RS(J,3) + FF(J,4)*RS(J,4) X 2 + FF(J,5)*RS(J,5) + FF(J,6)*RS(J,6) + FF(J,7)*RS(J,7) + FF X 3 (J,8)*RS(J,8) + FF(J,9)*RS(J,9) X 3400 CONTINUE X X DO 3600 J = I + 1, NMOL X IF (KC(J) .EQ. 0) THEN X RS(J,10) = XL(J,10)*XL(J,10) + YL(J,10)*YL(J,10) + ZL(J,10 X 1 )*ZL(J,10) X RL(J,10) = SQRT(RS(J,10)) X FF(J,10) = FUNEXP(AB1,B1,RL(J,10)) X X FTEMP(J,1) = FUNEXP(AB2,B2,RL(J,6)) X FF(J,6) = FF(J,6) + FTEMP(J,1) X RS(J,11) = XL(J,11)*XL(J,11) + YL(J,11)*YL(J,11) + ZL(J,11 X 1 )*ZL(J,11) X RL(J,11) = SQRT(RS(J,11)) X FF(J,11)=FUNEXP(AB3,B3,RL(J,11))-FUNEXP(AB4,B4,RL(J,11)) X X FTEMP(J,2) = FUNEXP(AB2,B2,RL(J,7)) X FF(J,7) = FF(J,7) + FTEMP(J,2) X RS(J,12) = XL(J,12)*XL(J,12) + YL(J,12)*YL(J,12) + ZL(J,12 X 1 )*ZL(J,12) X RL(J,12) = SQRT(RS(J,12)) X FF(J,12)=FUNEXP(AB3,B3,RL(J,12))-FUNEXP(AB4,B4,RL(J,12)) X X FTEMP(J,3) = FUNEXP(AB2,B2,RL(J,8)) X FF(J,8) = FF(J,8) + FTEMP(J,3) X RS(J,13) = XL(J,13)*XL(J,13) + YL(J,13)*YL(J,13) + ZL(J,13 X 1 )*ZL(J,13) X RL(J,13) = SQRT(RS(J,13)) X FF(J,13)=FUNEXP(AB3,B3,RL(J,13))-FUNEXP(AB4,B4,RL(J,13)) X X FTEMP(J,4) = FUNEXP(AB2,B2,RL(J,9)) X FF(J,9) = FF(J,9) + FTEMP(J,4) X RS(J,14) = XL(J,14)*XL(J,14) + YL(J,14)*YL(J,14) + ZL(J,14 X 1 )*ZL(J,14) X RL(J,14) = SQRT(RS(J,14)) X FF(J,14)=FUNEXP(AB3,B3,RL(J,14))-FUNEXP(AB4,B4,RL(J,14)) X ENDIF X 3600 CONTINUE X X DO 3800 J = I + 1, NMOL X IF (KC(J) .EQ. 0) VIRCPU(IDCPU) = VIRCPU(IDCPU) + FF(J,10)*RS( X 1 J,10) + FTEMP(J,1)*RS(J,6) + FF(J,11)*RS(J,11) + FTEMP(J,2) X 2 *RS(J,7) + FF(J,12)*RS(J,12) + FTEMP(J,3)*RS(J,8) + FF(J,13 X 3 )*RS(J,13) + FTEMP(J,4)*RS(J,9) + FF(J,14)*RS(J,14) X 3800 CONTINUE X X DO 4200 K = 1, 14 X DO 4000 JS = 1, JNUM X J = JSAVE(JS) X GG(JS,K) = FF(J,K)*XL(J,K) X 4000 CONTINUE X 4200 CONTINUE X X DO 4400 JS = 1, JNUM X J = JSAVE(JS) X JW1 = IW10 + (J-1)*NATOMS X JWO = IWO0 + (J-1)*NATOMS X JW2 = IW20 + (J-1)*NATOMS X G110 = GG(JS,10) + GG(JS,1)*C1 X G23 = GG(JS,2) + GG(JS,3) X G45 = GG(JS,4) + GG(JS,5) X FTEMP(JS,1) = G110 + GG(JS,11) + GG(JS,12) + C1*G23 X TT1 = GG(JS,1)*C2 X TT = G23*C2 + TT1 X FTEMP(JS,2) = GG(JS,6) + GG(JS,7) + GG(JS,13) + TT + GG(JS,4) X FTEMP(JS,3) = GG(JS,8) + GG(JS,9) + GG(JS,14) + TT + GG(JS,5) X TT = G45*C2 + TT1 X X X FXT(JWO,IDCPU)=FXT(JWO,IDCPU)-G110-GG(JS,13)-GG(JS,14)-C1*G45 X FXT(JW1,IDCPU) = FXT(JW1,IDCPU) - GG(JS,6) - GG(JS,8) - GG(JS, X 1 11) - TT - GG(JS,2) X FXT(JW2,IDCPU) = FXT(JW2,IDCPU) - GG(JS,7) - GG(JS,9) - GG(JS, X 1 12) - TT - GG(JS,3) X 4400 CONTINUE X X S1 = 0. X S2 = 0. X S3 = 0. X DO 4600 JS = 1, JNUM X S1 = S1 + FTEMP(JS,1) X S2 = S2 + FTEMP(JS,2) X S3 = S3 + FTEMP(JS,3) X 4600 CONTINUE X FX(IWO) = FX(IWO) + S1 X FX(IW1) = FX(IW1) + S2 X FX(IW2) = FX(IW2) + S3 X DO 5000 K = 1, 14 X DO 4800 JS = 1, JNUM X J = JSAVE(JS) X GG(JS,K) = FF(J,K)*YL(J,K) X 4800 CONTINUE X 5000 CONTINUE X X DO 5200 JS = 1, JNUM X J = JSAVE(JS) X JW1 = IW10 + (J-1)*NATOMS X JWO = IWO0 + (J-1)*NATOMS X JW2 = IW20 + (J-1)*NATOMS X X G110 = GG(JS,10) + GG(JS,1)*C1 X G23 = GG(JS,2) + GG(JS,3) X G45 = GG(JS,4) + GG(JS,5) X FTEMP(JS,1) = G110 + GG(JS,11) + GG(JS,12) + C1*G23 X TT1 = GG(JS,1)*C2 X TT = G23*C2 + TT1 X FTEMP(JS,2) = GG(JS,6) + GG(JS,7) + GG(JS,13) + TT + GG(JS,4) X FTEMP(JS,3) = GG(JS,8) + GG(JS,9) + GG(JS,14) + TT + GG(JS,5) X TT = G45*C2 + TT1 X FYT(JWO,IDCPU)=FYT(JWO,IDCPU)-G110-GG(JS,13)-GG(JS,14)-C1*G45 X FYT(JW1,IDCPU) = FYT(JW1,IDCPU) - GG(JS,6) - GG(JS,8) - GG(JS, X 1 11) - TT - GG(JS,2) X FYT(JW2,IDCPU) = FYT(JW2,IDCPU) - GG(JS,7) - GG(JS,9) - GG(JS, X 1 12) - TT - GG(JS,3) X 5200 CONTINUE X X S1 = 0. X S2 = 0. X S3 = 0. X DO 5400 JS = 1, JNUM X S1 = S1 + FTEMP(JS,1) X S2 = S2 + FTEMP(JS,2) X S3 = S3 + FTEMP(JS,3) X 5400 CONTINUE X FY(IWO) = FY(IWO) + S1 X FY(IW1) = FY(IW1) + S2 X FY(IW2) = FY(IW2) + S3 X X DO 5800 K = 1, 14 X DO 5600 JS = 1, JNUM X J = JSAVE(JS) X GG(JS,K) = FF(J,K)*ZL(J,K) X 5600 CONTINUE X 5800 CONTINUE X X DO 6000 JS = 1, JNUM X J = JSAVE(JS) X JW1 = IW10 + (J-1)*NATOMS X JWO = IWO0 + (J-1)*NATOMS X JW2 = IW20 + (J-1)*NATOMS X X G110 = GG(JS,10) + GG(JS,1)*C1 X G23 = GG(JS,2) + GG(JS,3) X G45 = GG(JS,4) + GG(JS,5) X FTEMP(JS,1) = G110 + GG(JS,11) + GG(JS,12) + C1*G23 X TT1 = GG(JS,1)*C2 X TT = G23*C2 + TT1 X FTEMP(JS,2) = GG(JS,6) + GG(JS,7) + GG(JS,13) + TT + GG(JS,4) X FTEMP(JS,3) = GG(JS,8) + GG(JS,9) + GG(JS,14) + TT + GG(JS,5) X TT = G45*C2 + TT1 X FZT(JWO,IDCPU)=FZT(JWO,IDCPU)-G110-GG(JS,13)-GG(JS,14)-C1*G45 X FZT(JW1,IDCPU) = FZT(JW1,IDCPU) - GG(JS,6) - GG(JS,8) - GG(JS, X 1 11) - TT - GG(JS,2) X FZT(JW2,IDCPU) = FZT(JW2,IDCPU) - GG(JS,7) - GG(JS,9) - GG(JS, X 1 12) - TT - GG(JS,3) X 6000 CONTINUE X X S1 = 0. X S2 = 0. X S3 = 0. X DO 6200 JS = 1, JNUM X S1 = S1 + FTEMP(JS,1) X S2 = S2 + FTEMP(JS,2) X S3 = S3 + FTEMP(JS,3) X 6200 CONTINUE X FZ(IWO) = FZ(IWO) + S1 X FZ(IW1) = FZ(IW1) + S2 X FZ(IW2) = FZ(IW2) + S3 X X 6400 CONTINUE X DO 6800 ID = 1, NCPUS X VIR = VIR + VIRCPU(ID) X DO 6600 I = 1, MAX X FX(I) = FX(I) + FXT(I,ID) X 6600 CONTINUE X 6800 CONTINUE X DO 7200 ID = 1, NCPUS X DO 7000 I = 1, MAX X FY(I) = FY(I) + FYT(I,ID) X 7000 CONTINUE X 7200 CONTINUE X DO 7600 ID = 1, NCPUS X DO 7400 I = 1, MAX X FZ(I) = FZ(I) + FZT(I,ID) X 7400 CONTINUE X 7600 CONTINUE X DO 7800 I = 1, NATMO, NATOMS X FX(I) = FX(I)*FHM X FY(I) = FY(I)*FHM X FZ(I) = FZ(I)*FHM X FX(I+2) = FX(I+2)*FHM X FY(I+2) = FY(I+2)*FHM X FZ(I+2) = FZ(I+2)*FHM X FX(I+1) = FX(I+1)*FOM X FY(I+1) = FY(I+1)*FOM X FZ(I+1) = FZ(I+1)*FOM X 7800 CONTINUE X X RETURN X END X SUBROUTINE SUB016 (DATA,EX) X X COMMON/DATA/NPTS,NSKIP,MTRN,MSKIP,ISIGN,LOG,IXSHFT C X COMPLEX DATA(NPTS), EXK, EXJ, H, EX(1), FACT X REAL HH(2) X EQUIVALENCE (HH,H) X X IERR = 0 X HH(1) = 0 X HH(2) = 0 X FACT = (1.2349E-5,1.8567E-5) X I2K = NPTS X IF (I2K .EQ. 1) GO TO 4000 X SGN1 = ISIGN X EXK = EX(1+IXSHFT) X IF (SGN1 .LT. 0.) EXK = CONJG(EXK) X 1000 CONTINUE X I2KP = I2K X I2K = I2K/2 X I2KS = I2K*NSKIP X IF (2*I2K .NE. I2KP) GO TO 4000 X JLI = I2K/2 + 1 X DO 3000 JL = 1, I2K X IF (JL - 1 .GT. 0) GO TO 1600 X EXJ = (1.,0.) X DO 1400 JJ = JL, NPTS, I2KP X DO 1200 MM = 1, MTRN X JS = (JJ-1)*NSKIP + (MM-1)*MSKIP + 1 X H = DATA(JS) - DATA(JS+I2KS) X DATA(JS) = (DATA(JS)+DATA(JS+I2KS))*FACT X DATA(JS+I2KS) = H*FACT X 1200 CONTINUE X 1400 CONTINUE X GO TO 2800 X 1600 CONTINUE X IF (JL - JLI .EQ. 0) GO TO 2200 C X EXJ = EXJ*EXK X DO 2000 JJ = JL, NPTS, I2KP X DO 1800 MM = 1, MTRN X JS = (JJ-1)*NSKIP + (MM-1)*MSKIP + 1 X H = DATA(JS) - DATA(JS+I2KS) X DATA(JS) = (DATA(JS)+DATA(JS+I2KS))*FACT X DATA(JS+I2KS) = (H*EXJ)*FACT X 1800 CONTINUE X 2000 CONTINUE X GO TO 2800 X 2200 CONTINUE X EXJ = CMPLX(0.,SGN1) X DO 2600 JJ = JL, NPTS, I2KP X DO 2400 MM = 1, MTRN X JS = (JJ-1)*NSKIP + (MM-1)*MSKIP + 1 X H = DATA(JS) - DATA(JS+I2KS) X DATA(JS) = (DATA(JS)+DATA(JS+I2KS))*FACT X DATA(JS+I2KS) = CMPLX((-SGN1*HH(2)),SGN1*HH(1))*FACT X 2400 CONTINUE X 2600 CONTINUE X 2800 CONTINUE X 3000 CONTINUE C X EXK = EXK*EXK X IF (I2K - 1 .GT. 0) GO TO 1000 X IF (NPTS .LE. 2) GO TO 4000 X NPTSD2 = NPTS/2 X JMIN = 0 X JMAX = NPTS - 4 X JREV = 0 X DO 3800 J = JMIN, JMAX, 2 X I2K = NPTSD2 X JREV2 = JREV + I2K X IF (JREV2 - (J+1) .LE. 0) GO TO 3600 X DO 3200 MM = 1, MTRN X JREV2S = (MM-1)*MSKIP + JREV2*NSKIP X JS = (MM-1)*MSKIP + (J+1)*NSKIP X H = DATA(JREV2S+1) X DATA(JREV2S+1) = DATA(JS+1)*FACT X DATA(JS+1) = H*FACT X 3200 CONTINUE X IF (JREV - J .LE. 0) GO TO 3600 X DO 3400 MM = 1, MTRN X JREVS = JREV*NSKIP + (MM-1)*MSKIP X JS = J*NSKIP + (MM-1)*MSKIP X H = DATA(JREVS+1) X DATA(JREVS+1) = DATA(JS+1)*FACT X DATA(JS+1) = H*FACT X 3400 CONTINUE X 3600 CONTINUE X 3800 CONTINUE C X 4000 CONTINUE C X NFTVMT = NFTVMT + 1 C X RETURN X END X SUBROUTINE SUB017 C X PARAMETER (NNX=44) X PARAMETER (NNY=44) X PARAMETER (NNZ=4) X PARAMETER (MXC=11) X PARAMETER (MXP=3) X PARAMETER (MXW=5) X PARAMETER (NXYN=NNX*NNY) X PARAMETER (NBL=NXYN*NNZ) X PARAMETER (NBLW=NBL+MXW) X PARAMETER (N=MXC) X PARAMETER (NY=NNY) X PARAMETER (NX=NNX) X PARAMETER (NZ=NNZ) X PARAMETER (NXM1=NX-1) X PARAMETER (NWELL=5) C ---------------------------------------------------------------- X COMMON /DATA/ CTOT(NBLW,MXC),C(NBLW,MXC,MXP) X * ,CSE(NBLW),S(NBLW,MXP) X COMMON /DATA/ BRK,CRK,VIS1,VIS2,AP1,AP2,AP3,GAMMAC,EPHI4,EPHI3, X * GAMHF,SSLOPE,POWN,CSE1,ALPHA1,ALPHA2,ALPHA3,ALPHA4,ALPHA5, X * BETAP X COMMON /DATA/ VIS(NBL,MXP),RPERM(NBLW,MXP) X *,PERMX(NBL),PERMY(NBL),PERMZ(NBL), X *QI(MXW,MXP),QT(MXW,NNZ),Q(MXW,NNZ,MXP),PWF(MXW) X COMMON /DATA/ EL(NBL),DX(NNX),DY(NNY), X * DZ(NNZ) X COMMON /DATA/ TRSX(NBL,MXP),TRSY(NBL,MXP), X *TRSZ(NBL,MXP),TX(NBL),TY(NBL), X *TZ(NBL),CONVX(NBL,MXC,MXP),CONVY(NBL,MXC,MXP), X *CONVZ(NBL,MXC,MXP),VX(NBL,MXP),VY(NBL,MXP), X *VZ(NBL,MXP) X COMMON /DATA/ POR(NBL),RKF(NBL,MXP) X COMMON /DATA/ RW(MXW),SWELL(MXW),PI(MXW,NNZ),TTM(NNZ), X $ TM(MXW,NNZ),PTM(MXW,NNZ,MXP),GAMMAT(MXW,NNZ) X COMMON /DATA/ C6JO(NBLW),C6ADSS(NBLW), X * C6HATS(NBLW),QV,XKC,XKS,EQW X COMMON /DATA/ CUMI(MXC),CUMP(MXC),OIP(MXC), X * T,TINJ,ICNT,WHPV,PRF X COMMON /DATA/FOREC,FOREC1,RESPV,RESATK(NNZ),BTO(MXW), X * VB(NBL) X COMMON /DATA/ TITLE(30),RELERR(MXC),SWI X COMMON /DATA/ D(MXC,MXP),ALPHAL(MXP),ALPHAT(MXP) X COMMON /DATA/ DCMAX,DISPC X LOGICAL LWKSP1 X COMMON /DATA/ LWKSP1(NBL) X COMMON /DATA/ COE1(NBLW),COE2(NBLW),COE3(NBLW), X * COE4(NBLW),COE5(NBLW) X COMMON /DATA/ COE6(NBLW),COE7(NBLW),COEX(NBLW),COEYU(NBLW), X * COEZT(NBLW),DDX(NBLW),DDY(NBLW),DDZ(NBLW), X * DC(NBLW) X COMMON /DATA/ TK(3),TTRP(3),TREC(3),TBT(MXW,3), X * RDC(3),TRD(3),RET(3) X COMMON /DATA/ AG1,AG2,CRG,AGK,BGK X COMMON /DATA/ A15D,B15D,C15DSS(NBL),C14DSS(NBL) X * ,QW(NBL),TC14DS(NBL) X COMMON /DATA/ AK1,AK2,SCR,GKIN(MXC),DT C X COMMON /INTS/ IPERM,IMES X COMMON /INTS/ NGP,NGP1,NGP2,NGP3,NGP4,NG,IREACT, X $ IW(MXW),JW(MXW),IJKW(MXW),KWELL(MXW,NNZ),ICF(MXC) C X DIMENSION OP(MXC),TEMP(MXC) C X IJK = 0 X DO 1400 K = 1, NZ X DO 1200 J = 1, NY X DO 1000 I = 1, NX X IJK = IJK + 1 X DDX(IJK) = DX(I) X DDY(IJK) = DY(J) X DDZ(IJK) = DZ(K) X 1000 CONTINUE X 1200 CONTINUE X 1400 CONTINUE X IF (IMES .EQ. 2) DCMAX = 0.0 X DO 20000 KC = 1, N X IF (ICF(KC) .EQ. 0) GO TO 19800 X EPHI1 = 1.0 X IF (KC .EQ. 3) EPHI1 = EPHI3 X IF (KC .EQ. 4) EPHI1 = EPHI4 X DO 1600 I = 1, NBL X COEX(I) = 0.0 X COEYU(I) = 0.0 X COEZT(I) = 0.0 X 1600 CONTINUE X DO 13200 L = 1, 3 X ALPHAX = ALPHAL(L) - ALPHAT(L) C X IF (NY .GT. 1) THEN C X DO 1800 I = NX + 1, NBL - NX X COE1(I) = 0.25*(VY(I,L)+VY(I+1,L)+VY(I-NX,L)+VY(I+1-NX,L X 1 )) X 1800 CONTINUE X DO 2400 K = 1, NZ X IBGN = (K-1)*NXYN + 1 X IEND = IBGN + NXM1 X DO 2000 I = IBGN, IEND X COE1(I) = 0.5*(VY(I,L)+VY(I+1,L)) X 2000 CONTINUE X IBGN = K*NXYN - NXM1 X IEND = IBGN + NXM1 X DO 2200 I = IBGN, IEND X COE1(I) = 0.5*(VY(I-NX,L)+VY(I-NX+1,L)) X 2200 CONTINUE X 2400 CONTINUE X ELSE X DO 2600 I = 1, NBL X COE1(I) = 0.0 X 2600 CONTINUE X ENDIF X IF (NZ .GT. 1) THEN X DO 2800 I = 1, NXYN X COE2(I) = 0.5*(VZ(I,L)+VZ(I+1,L)) X 2800 CONTINUE X IBGN = (NZ-1)*NXYN + 1 X DO 3000 I = IBGN, NBL X COE2(I) = 0.5*(VZ(I-NXYN,L)+VZ(I-NXYN+1,L)) X 3000 CONTINUE X IF (NZ .GT. 2) THEN X DO 3200 I = NXYN + 1, NBL - NXYN X COE2(I) = 0.25*(VZ(I,L)+VZ(I+1,L)+VZ(I-NXYN,L)+VZ(I- X 1 NXYN+1,L)) X 3200 CONTINUE X ENDIF X ELSE X DO 3400 I = 1, NBL X COE2(I) = 0.0 X 3400 CONTINUE X ENDIF X DO 3600 I = 1, NBL X COE3(I) = VX(I,L)**2 X COE4(I) = COE1(I)**2 X COE5(I) = COE2(I)**2 X COE6(I) = SQRT(COE3(I)+COE4(I)+COE5(I)) X IF (COE6(I) .LE. .0001) COE6(I) = 1.E199 X X COE6(I) = 1.0/COE6(I) X 3600 CONTINUE X DO 3800 I = 1, NBL - 1 X COE7(I) = 0.5*(S(I,L)*POR(I)+S(I+1,L)*POR(I+1)) X 3800 CONTINUE C X DO 4000 I = 1, NBL X COE7(I) = COE7(I)*D(KC,L) + (ALPHAL(L)*COE3(I)+ALPHAT(L)*( X 1 COE4(I)+COE5(I)))*COE6(I) X COE1(I) = ALPHAX*ABS(COE1(I)*VX(I,L))*COE6(I) X COE2(I) = ALPHAX*ABS(COE2(I)*VX(I,L))*COE6(I) X 4000 CONTINUE X DO 4200 I = 1, NBL - 1 X COE3(I) = 2.*(.5*DISPC*DDX(I)*ABS(VX(I,L))-COE7(I))*(C(I+1, X 1 KC,L)-C(I,KC,L))/(DDX(I+1)+DDX(I)) X 4200 CONTINUE X IF (NY .GT. 2) THEN X DO 4400 I = NX + 1, NBL - NX - 1 X COE4(I) = 0.25*COE1(I)*(C(I+NX,KC,L)+C(I+NX+1,KC,L)-C(I- X 1 NX,KC,L)-C(I-NX+1,KC,L))/DDY(I) X 4400 CONTINUE X DO 5000 K = 1, NZ X IBGN = (K-1)*NXYN + 1 X IEND = IBGN + NXM1 X DO 4600 I = IBGN, IEND X COE4(I) = 0.0 X 4600 CONTINUE X IBGN = K*NXYN - NXM1 X IEND = IBGN + NXM1 X DO 4800 I = IBGN, IEND X COE4(I) = 0.0 X 4800 CONTINUE X 5000 CONTINUE X ELSE X DO 5200 I = 1, NBL X COE4(I) = 0.0 X 5200 CONTINUE X ENDIF X DO 5400 I = 1, NBL X COE5(I) = 0.0 X 5400 CONTINUE X IF (NZ .GT. 2) THEN X DO 5600 I = NXYN + 1, NBL - NXYN - 1 X COE5(I) = 0.25*COE2(I)*(C(I+NXYN,KC,L)+C(I+NXYN+1,KC,L)- X 1 C(I-NXYN,KC,L)-C(I-NXYN+1,KC,L))/DDZ(I) X 5600 CONTINUE X ENDIF C X DO 5800 I = 1, NBL - 1 X COEX(I+1) = COEX(I+1) + COE3(I) - COE4(I) - COE5(I) X 5800 CONTINUE C X IF (NY .GT. 1) THEN C X DO 6000 I = 2, NBL - NX X COE1(I) = 0.25*(VX(I,L)+VX(I+NX,L)+VX(I-1,L)+VX(I+NXM1,L X 1 )) X 6000 CONTINUE X DO 6200 I = 1, NBL - NX, NX X COE1(I) = 0.5*(VX(I,L)+VX(I+NX,L)) X 6200 CONTINUE X DO 6400 I = NX, NBL - NX, NX X COE1(I) = 0.5*(VX(I-1,L)+VX(I-1+NX,L)) X 6400 CONTINUE X IF (NZ .GT. 1) THEN X DO 6600 I = 1, NXYN X COE2(I) = 0.5*(VZ(I,L)+VZ(I+NX,L)) X 6600 CONTINUE X DO 6800 I = (NZ-1)*NXYN + 1, NBL X COE2(I) = 0.5*(VZ(I-NXYN,L)+VZ(I-NXYN+NX,L)) X 6800 CONTINUE X IF (NZ .GT. 2) THEN X DO 7000 I = NXYN + 1, (NZ-1)*NXYN X COE2(I) = 0.25*(VZ(I,L)+VZ(I+NX,L)+VZ(I-NXYN,L)+ X 1 VZ(I-NXYN+NX,L)) X 7000 CONTINUE X ENDIF X ELSE X DO 7200 I = 1, NBL X COE2(I) = 0.0 X 7200 CONTINUE X ENDIF X DO 7400 I = 1, NBL X COE3(I) = VY(I,L)**2 X COE4(I) = COE1(I)**2 X COE5(I) = COE2(I)**2 X COE6(I) = SQRT(COE3(I)+COE4(I)+COE5(I)) X IF (COE6(I) .LE. .0001) COE6(I) = 1.E199 X COE6(I) = 1.0/COE6(I) X 7400 CONTINUE X DO 7600 I = 1, NBL - NX X COE7(I) = 0.5*(S(I,L)*POR(I)+S(I+NX,L)*POR(I+NX)) X 7600 CONTINUE X DO 7800 I = 1, NBL X COE7(I) = COE7(I)*D(KC,L) + (ALPHAL(L)*COE3(I)+ALPHAT(L) X 1 *(COE4(I)+COE5(I)))*COE6(I) X COE1(I) = ALPHAX*ABS(COE1(I)*VY(I,L))*COE6(I) X COE2(I) = ALPHAX*ABS(COE2(I)*VY(I,L))*COE6(I) X 7800 CONTINUE C X DO 8000 I = 1, NBL - NX X COE3(I) = 2.*(.5*DISPC*DDY(I)*ABS(VY(I,L))-COE7(I))*(C(I X 1 +NX,KC,L)-C(I,KC,L))/(DDY(I+NX)+DDY(I)) X 8000 CONTINUE X DO 8200 I = 2, NBL - NX - 1 X COE4(I) = 0.25*COE1(I)*(C(I+1,KC,L)+C(I+1+NX,KC,L)-C(I-1 X 1 ,KC,L)-C(I-1+NX,KC,L))/DDX(I) X 8200 CONTINUE X DO 8400 I = 1, NBL, NX X COE4(I) = 0.0 X 8400 CONTINUE X DO 8600 I = NX, NBL, NX X COE4(I) = 0.0 X 8600 CONTINUE X DO 8800 I = 1, NBL X COE5(I) = 0.0 X 8800 CONTINUE X IF (NZ .GT. 2) THEN X DO 9000 I = NXYN + 1, NBL - NXYN - NX X COE5(I) = 0.25*COE2(I)*(C(I+NXYN,KC,L)+C(I+NXYN+NX, X 1 KC,L)-C(I-NXYN,KC,L)-C(I-NXYN+NX,KC,L))/DDZ(I) X 9000 CONTINUE X ENDIF C X DO 9200 I = 1, NBL - NX X COEYU(I+NX) = COEYU(I+NX) + COE3(I) - COE4(I) - COE5(I) X 9200 CONTINUE X ENDIF C X IF (NZ .GT. 1) THEN C X DO 9400 I = 2, NBL - NXYN X COE1(I) = 0.25*(VX(I,L)+VX(I-1,L)+VX(I+NXYN,L)+VX(I+NXYN X 1 -1,L)) X 9400 CONTINUE X DO 9600 I = 1, NBL - NXYN, NX X COE1(I) = 0.5*(VX(I,L)+VX(I+NXYN,L)) X 9600 CONTINUE X DO 9800 I = NX, NBL - NXYN, NX X COE1(I) = 0.5*(VX(I-1,L)+VX(I+NXYN-1,L)) X 9800 CONTINUE X DO 10000 I = NX + 1, NBL - NXYN - NX X COE2(I) = 0.25*(VY(I,L)+VY(I-NX,L)+VY(I+NXYN,L)+VY(I+ X 1 NXYN-NX,L)) X10000 CONTINUE X DO 10600 K = 1, NZ - 1 X IBGN = (K-1)*NXYN + 1 X IEND = IBGN + NXM1 X DO 10200 I = IBGN, IEND X COE2(I) = 0.5*(VY(I,L)+VY(I+NXYN,L)) X10200 CONTINUE X IBGN = K*NXYN - NXM1 X IEND = IBGN + NXM1 X DO 10400 I = IBGN, IEND X COE2(I) = 0.5*(VY(I-NX,L)+VY(I+NXYN-NX,L)) X10400 CONTINUE X10600 CONTINUE X DO 10800 I = 1, NBL - NXYN X COE3(I) = VZ(I,L)**2 X COE4(I) = COE1(I)**2 X COE5(I) = COE2(I)**2 X COE6(I) = SQRT(COE3(I)+COE4(I)+COE5(I)) X IF (COE6(I) .LE. .0001) COE6(I) = 1.E199 X COE6(I) = 1.0/COE6(I) X10800 CONTINUE X DO 11000 I = 1, NBL - NXYN X COE7(I) = 0.5*(S(I,L)*POR(I)+S(I+NXYN,L)*POR(I+NXYN)) X11000 CONTINUE X DO 11200 I = 1, NBL - NXYN X COE7(I) = COE7(I)*D(KC,L) + (ALPHAL(L)*COE3(I)+ALPHAT(L) X 1 *(COE4(I)+COE5(I)))*COE6(I) X COE1(I) = ALPHAX*ABS(COE1(I)*VZ(I,L))*COE6(I) X COE2(I) = ALPHAX*ABS(COE2(I)*VZ(I,L))*COE6(I) X11200 CONTINUE C X DO 11400 I = 1, NBL - NXYN X COE3(I) = 2.*(.5*DISPC*DDZ(I)*ABS(VZ(I,L))-COE7(I))*(C(I X 1 +NXYN,KC,L)-C(I,KC,L))/(DDZ(I+NXYN)+DDZ(I)) X11400 CONTINUE X DO 11600 I = 2, NBL - NXYN - 1 X COE4(I) = 0.25*COE1(I)*(C(I+1,KC,L)+C(I+NXYN+1,KC,L)-C(I X 1 -1,KC,L)-C(I+NXYN-1,KC,L))/DDX(I) X11600 CONTINUE X DO 11800 I = 1, NBL - NXYN, NX X COE4(I) = 0.0 X11800 CONTINUE X DO 12000 I = NX, NBL - NXYN, NX X COE4(I) = 0.0 X12000 CONTINUE X DO 12200 I = NX + 1, NBL - NXYN - NX X COE5(I) = 0.25*COE2(I)*(C(I+NX,KC,L)+C(I+NX+NXYN,KC,L)-C X 1 (I-NX,KC,L)-C(I+NXYN-NX,KC,L))/DDY(I) X12200 CONTINUE X DO 12800 K = 1, NZ - 1 X IBGN = (K-1)*NXYN + 1 X IEND = IBGN + NXM1 X DO 12400 I = IBGN, IEND X COE5(I) = 0.0 X12400 CONTINUE X IBGN = K*NXYN - NXM1 X IEND = IBGN + NXM1 X DO 12600 I = IBGN, IEND X COE5(I) = 0.0 X12600 CONTINUE X12800 CONTINUE C X DO 13000 I = 1, NBL - NXYN X COEZT(I+NXYN)=COEZT(I+NXYN)+COE3(I)-COE4(I)-COE5(I) X13000 CONTINUE X ENDIF X13200 CONTINUE X DO 13400 I = 1, NBL, NX X COEX(I) = 0.0 X13400 CONTINUE X DO 13800 K = 1, NZ X IBGN = (K-1)*NXYN + 1 X IEND = IBGN + NXM1 X DO 13600 I = IBGN, IEND X COEYU(I) = 0.0 X13600 CONTINUE X13800 CONTINUE X DO 14000 I = 1, NXYN X COEZT(I) = 0.0 X14000 CONTINUE X DO 14200 I = 1, NBL X CONVX(I,KC,1) = CONVX(I,KC,1) + CONVX(I,KC,2) + CONVX(I,KC,3) X CONVY(I,KC,1) = CONVY(I,KC,1) + CONVY(I,KC,2) + CONVY(I,KC,3) X CONVZ(I,KC,1) = CONVZ(I,KC,1) + CONVZ(I,KC,2) + CONVZ(I,KC,3) X14200 CONTINUE X DO 14400 I = 1, NBL X COE4(I) = 0.0 X COE5(I) = 0.0 X14400 CONTINUE X COE3(1) = (COEX(2)+CONVX(1,KC,1))/DDX(1) X DO 14600 I = 2, NBL - 1 X COE3(I) = (COEX(I+1)-COEX(I)+CONVX(I,KC,1)-CONVX(I-1,KC,1))/ X 1 DDX(I) X14600 CONTINUE X COE3(NBL) = ((-COEX(NBL))+CONVX(NBL,KC,1)-CONVX(NBL-1,KC,1))/DDX X 1 (NBL) X IF (NY .GT. 1) THEN X DO 14800 I = 1, NX X COE4(I) = (COEYU(I+NX)+CONVY(I,KC,1))/DDY(I) X14800 CONTINUE X DO 15000 I = NX + 1, NBL - NX X COE4(I) = (COEYU(I+NX)-COEYU(I)+CONVY(I,KC,1)-CONVY(I-NX, X 1 KC,1))/DDY(I) X15000 CONTINUE X DO 15200 I = NBL - NX + 1, NBL X COE4(I) = ((-COEYU(I))+CONVY(I,KC,1)-CONVY(I-NX,KC,1))/DDY X 1 (I) X15200 CONTINUE X ENDIF X IF (NZ .GT. 1) THEN X IF (NZ .EQ. 2) THEN X DO 15400 I = 1, NXYN X COE5(I) = (COEZT(I+NXYN)+CONVZ(I,KC,1))/DDZ(I) X15400 CONTINUE X DO 15600 I = NXYN + 1, NBL X COE5(I) = ((-COEZT(I))+CONVZ(I,KC,1)-CONVZ(I-NXYN,KC,1 X 1 ))/DDZ(I) X15600 CONTINUE X ELSE X DO 15800 I = NXYN + 1, NBL - NXYN X COE5(I) = (COEZT(I+NXYN)-COEZT(I)+CONVZ(I,KC,1)-CONVZ( X 1 I-NXYN,KC,1))/DDZ(I) X15800 CONTINUE X DO 16000 I = NBL - NXYN + 1, NBL X COE5(I) = ((-COEZT(I))+CONVZ(I,KC,1)-CONVZ(I-NXYN,KC,1 X 1 ))/DDZ(I) X16000 CONTINUE X ENDIF X ENDIF C X IF (KC .GT. 8) THEN X IT = KC - 8 X IF (RET(IT) .NE. 0.) THEN X DO 16200 I = 1, NBL X CON = 1./(1.+RET(IT)/S(I,1)) X COE3(I) = COE3(I)*CON X COE4(I) = COE4(I)*CON X COE5(I) = COE5(I)*CON X16200 CONTINUE X ENDIF X ENDIF X DO 16400 I = 1, NBL X DC(I) = (-DT/EPHI1)*(COE3(I)+COE4(I)+COE5(I))/POR(I) X16400 CONTINUE C X DO 17000 M = 1, NWELL X I = IW(M) X J = JW(M) X DO 16800 K = 1, NZ X IJK1 = (K-1)*NXYN + (J-1)*NX + I X DO 16600 L = 1, 3 X IF (Q(M,K,L) .GT. 0.0) THEN X CW = C(IJKW(M),KC,L) X CUMI(KC) = CUMI(KC) + Q(M,K,L)*CW*DT X ELSE X CW = C(IJK1,KC,L) X IF (NY .NE. 1) THEN X TERM0 = 28.*CW X TERM1 = CW X IF (J .NE. 1) TERM1 = C(IJK1-NX,KC,L) X TERM2 = CW X IF (I .NE. 1) TERM2 = C(IJK1-1,KC,L) X TERM3 = CW X IF (I .NE. NX) TERM3 = C(IJK1+1,KC,L) X TERM4 = CW X IF (J .NE. NY) TERM4 = C(IJK1+NX,KC,L) X CW = (TERM0-TERM1-TERM2-TERM3-TERM4)/24. X ENDIF X CUMP(KC) = CUMP(KC) + Q(M,K,L)*CW*DT X ENDIF X DC(IJK1) = DC(IJK1) + (DT*CW/EPHI1)*Q(M,K,L)/VB(IJK1) X16600 CONTINUE X16800 CONTINUE X17000 CONTINUE C X IF (KC .GT. 8) THEN X KCT = KC - 8 X CONS = RDC(KCT)*DT X IF (CONS .EQ. 0.0) GO TO 18000 X DO 17200 I = 1, NBL X COE1(I) = 0.0 X17200 CONTINUE X DO 17600 L = 1, 3 X DO 17400 I = 1, NBL X COE1(I) = COE1(I) + CONS*C(I,KC,L) X17400 CONTINUE X17600 CONTINUE X TR = 0.0 X DO 17800 I = 1, NBL X DC(I) = DC(I) - COE1(I) X TR = TR + COE1(I)*VB(I) X17800 CONTINUE X TRD(KCT) = TRD(KCT) + TR X18000 CONTINUE X ENDIF C X IF (IREACT.EQ.1 .AND. NG.NE.0) THEN X IF (KC .EQ. 4) THEN X DO 18200 I = 1, NBL X COE1(I) = DT*(1./(SCR+1.))*AK2*1.E12*(C(I,NGP3,1)* X 1 1.E-4*C(I,4,1))**2 X DC(I) = DC(I) - COE1(I) X GKIN(4) = GKIN(4) + COE1(I)*VB(I) X18200 CONTINUE X ELSE IF (KC .EQ. NGP1) THEN X DO 18400 I = 1, NBL X COE1(I) = DT*AK1*C(I,NGP1,1)*C(I,NGP2,1) X DC(I) = DC(I) - COE1(I) X GKIN(NGP1) = GKIN(NGP1) + COE1(I)*VB(I) X18400 CONTINUE X ELSE IF (KC .EQ. NGP2) THEN X DO 18600 I = 1, NBL X COE1(I) = DT*AK1*2.111*C(I,NGP1,1)*C(I,NGP2,1) X DC(I) = DC(I) - COE1(I) X GKIN(NGP2) = GKIN(NGP2) + COE1(I)*VB(I) X18600 CONTINUE X ELSE IF (KC .EQ. NGP3) THEN X DO 18800 I = 1, NBL X COE1(I) = DT*(0.48*AK1*C(I,NGP1,1)*C(I,NGP2,1)-AK2*( X 1 SCR/(SCR+1.))*(C(I,4,1)*C(I,NGP3,1)*1.E4)**2) X DC(I) = DC(I) + COE1(I) X GKIN(NGP3) = GKIN(NGP3) + COE1(I)*VB(I) X18800 CONTINUE X ELSE IF (KC .EQ. NGP4) THEN X DO 19000 I = 1, NBL X COE1(I) = DT*AK2*(C(I,NGP3,1)*C(I,4,1)*1.E4)**2 X DC(I) = DC(I) + COE1(I) X GKIN(NGP4) = GKIN(NGP4) + COE1(I)*VB(I) X19000 CONTINUE X ENDIF X ENDIF X DO 19200 I = 1, NBL X CTOT(I,KC) = CTOT(I,KC) + DC(I) X19200 CONTINUE X IF (KC.LE.3 .AND. IMES.EQ.2) THEN X JDC = ISAMAX(NBL,DC,1) X DCMAX = AMAX1(DCMAX,ABS(DC(JDC))) X ENDIF X DO 19400 I = 1, NBL X LWKSP1(I) = CTOT(I,KC) .LE. (-1.E-2) X19400 CONTINUE C X SUM = 0.0 X DO 19600 I = 1, NBL X SUM = SUM + CTOT(I,KC)*VB(I) X19600 CONTINUE X OP(KC) = SUM*EPHI1 X IF (KC .GT. 8) OP(KC) = OP(KC) + TRD(KC-8)*EPHI1 X19800 CONTINUE X20000 CONTINUE X SUM = 0.0 X DO 20200 I = 1, NBL X SUM = SUM + (C6HATS(I)+C6ADSS(I))*VB(I) X20200 CONTINUE X OP(6) = OP(6) + SUM*EPHI1 X DO 20400 KC = 1, N X RELERR(KC) = 0.0 X TEMP(KC) = CUMI(KC) + OIP(KC) X IF (TEMP(KC) .GT. 1.E-7) RELERR(KC) = ABS(TEMP(KC)+CUMP(KC)-OP( X 1 KC))/TEMP(KC) X20400 CONTINUE X RETURN X END X SUBROUTINE SUB018(A,U,F,Z) X PARAMETER (NXY=100000,NX=100000) X DIMENSION A(NXY,7),U(NXY),F(NXY),Z(NXY) X Z(1) = 0.0 X KE = NXY - NX + 2 X NU = NX - 2 X X DO 1200 KK = 2, KE, 65535 X KKE = (KK-1) + MIN0(65535,KE-(KK-1)) X DO 1000 K = KK, KKE X Z(K) = A(K,3)*A(K-1,6)*U(K+NU) X 1000 CONTINUE X 1200 CONTINUE X X KB = KE + 1 X DO 1400 K = KB, NXY X Z(K) = 0.0 X 1400 CONTINUE X NA = (-NX) + 1 X NU = NA + 1 X X DO 1600 K = 1, NX X Z(K) = Z(K) + F(K) X 1600 CONTINUE X X DO 2000 KK = NX + 1, NXY, 65535 X KKE = (KK-1) + MIN0(65535,NXY-(KK-1)) X DO 1800 K = KK, KKE X Z(K) = (Z(K)+F(K)) + A(K,2)*A(K+NA,5)*U(K+NU) X 1800 CONTINUE X 2000 CONTINUE X X RETURN X END X SUBROUTINE SUB019 X PARAMETER (MR=100,MT=16,MX=16) X PARAMETER (PI=3.141592653589793) X COMMON /DATA/ X 1 VR(MR+1,MT+1,6*MX),VC(MR+1,MT+1,6*MX),WR(MR+1,MT+1,6*MX), X 1 WC(MR+1,MT+1,6*MX) X COMPLEX CVEC(MR+1),CFAC X DO 1400 J = 1, MT + 1 X DO 1200 K = 1, 2*MX X X J1 = J - 1 X K1 = K - 1 X IF (K1 .GT. MX) K1 = K1 - 2*MX X CFAC = EXP((0.,1.)*(J1*PI/(2*MT)+K1*PI/(2*MX))) X DO 1000 I = 1, MR + 1 X X CVEC(I) = VR(I,J,K) + (0.,1.)*VC(I,J,K) X CVEC(I) = CFAC*CVEC(I) X VR(I,J,K) = CVEC(I) X VC(I,J,K) = (0.,-1.)*CVEC(I) X X CVEC(I) = VR(I,J,K+2*MX) + (0.,1.)*VC(I,J,K+2*MX) X CVEC(I) = CFAC*CVEC(I) X VR(I,J,K+2*MX) = CVEC(I) X VC(I,J,K+2*MX) = (0.,-1.)*CVEC(I) X X CVEC(I) = VR(I,J,K+4*MX) + (0.,1.)*VC(I,J,K+4*MX) X CVEC(I) = CFAC*CVEC(I) X VR(I,J,K+4*MX) = CVEC(I) X VC(I,J,K+4*MX) = (0.,-1.)*CVEC(I) X 1000 CONTINUE X 1200 CONTINUE X 1400 CONTINUE X X DO 2000 J = 1, MT + 1 X DO 1800 K = 1, 2*MX X J1 = J - 1 X K1 = K - 1 X IF (K1 .GT. MX) K1 = K1 - 2*MX X X CFAC = EXP((0.,-1.)*(J1*PI/(2*MT)+K1*PI/(2*MX))) X DO 1600 I = 1, MR + 1 X CVEC(I) = WR(I,J,K) + (0.,1.)*WC(I,J,K) X CVEC(I) = CFAC*CVEC(I) X WR(I,J,K) = CVEC(I) X WC(I,J,K) = (0.,-1.)*CVEC(I) X X CVEC(I) = WR(I,J,K+2*MX) + (0.,1.)*WC(I,J,K+2*MX) X CVEC(I) = CFAC*CVEC(I) X WR(I,J,K+2*MX) = CVEC(I) X WC(I,J,K+2*MX) = (0.,-1.)*CVEC(I) X X CVEC(I) = WR(I,J,K+4*MX) + (0.,1.)*WC(I,J,K+4*MX) X CVEC(I) = CFAC*CVEC(I) X WR(I,J,K+4*MX) = CVEC(I) X WC(I,J,K+4*MX) = (0.,-1.)*CVEC(I) X X 1600 CONTINUE X 1800 CONTINUE X 2000 CONTINUE X RETURN X END X SUBROUTINE SUB020(S,KKP,KKR) X PARAMETER(JKLMAX=40,IMX=29791) X COMMON/DATA/Q(5,IMX), X(IMX), Y(IMX), Z(IMX) X COMMON/DATA/R(JKLMAX), RTXYZ(4,JKLMAX), DJ(IMX) X DIMENSION KKP(JKLMAX),KKR(JKLMAX),S(5,IMX) X DATA JMAX/JKLMAX/,KMAX/JKLMAX/,LMAX/JKLMAX/ X DATA JM/JKLMAX/,KM/JKLMAX/,LM/JKLMAX/ X DATA JMM/JKLMAX/,KMM/JKLMAX/LMM/JKLMAX/ X DATA KEND2/2/,KENDM/JKLMAX/,KK/0/,LL/0/ X DATA SMU/1.0/,SMUIM/1.0/ X DO 2000 L = 2, LM X DO 1800 K = KEND2, KENDM X IKL = (K-1)*KK + (L-1)*LL X DO 1200 N = 1, 4 X DO 1000 J = 3, JMM X I = IKL + J X SMJ = SMU*DJ(I) X D4Q=Q(N,I+2)-4.*Q(N,I+1)+6.*Q(N,I)-4.*Q(N,I-1)+Q(N,I-2) X S(N,I) = S(N,I) - SMJ*D4Q X 1000 CONTINUE X 1200 CONTINUE C. J=2 X I = IKL + 2 X SMJ = SMU*DJ(I) X DO 1400 N = 1, 4 X D4Q=Q(N,I-1)-4.*Q(N,I)+6.*Q(N,I+1)-4.*Q(N,I+2)+Q(N,I+3) X S(N,I) = S(N,I) - SMJ*D4Q X 1400 CONTINUE X I = IKL + JM X SMJ = SMU*DJ(I) X DO 1600 N = 1, 4 X D2Q = Q(N,I+1) - 2.*Q(N,I) + Q(N,I-1) X S(N,I) = S(N,I) + SMJ*D2Q X 1600 CONTINUE X 1800 CONTINUE X 2000 CONTINUE X 2200 CONTINUE X DO 5600 L = 2, LM X DO 5400 J = 2, JM X IJL = J + (L-1)*LL X DO 2600 N = 1, 4 X K = 2 X I = IJL + (K-1)*KK X DO 2400 K = 3, KMM X I = I + KK X SMJ = SMU*DJ(I) X D4Q = Q(N,I-KK-KK) - 4.*(Q(N,I-KK)+Q(N,I+KK)) + 6.*Q(N,I) X 1 + Q(N,I+KK+KK) X S(N,I) = S(N,I) - SMJ*D4Q X 2400 CONTINUE X 2600 CONTINUE X I = IJL + KK X SMJ = SMU*DJ(I) X DO 2800 N = 1, 4 X D4Q = Q(N,I-KK) - 4.*Q(N,I) + 6.*Q(N,I+KK) - 4.*Q(N,I+KK+KK) X 1 + Q(N,I+KK+KK+KK) X S(N,I) = S(N,I) - SMJ*D4Q X 2800 CONTINUE X I = IJL + (KM-1)*KK X SMJ = SMU*DJ(I) X DO 3000 N = 1, 4 X D4Q = Q(N,I+KK) - 4.*Q(N,I) + 6.*Q(N,I-KK) - 4.*Q(N,I-KK-KK) X 1 + Q(N,I-KK-KK-KK) X S(N,I) = S(N,I) - SMJ*D4Q X 3000 CONTINUE X GO TO 5200 X 3200 CONTINUE C. K=2 X I = IJL + KK X SMJ = SMU*DJ(I) X DO 3400 N = 1, 4 X D4Q = Q(N,I+KK+KK) - 4.*(Q(N,I+KK)+Q(N,I-KK)) + 7.*Q(N,I) X IF (N .EQ. 3) D4Q = D4Q - 2.*Q(3,I) X S(N,I) = S(N,I) - SMJ*D4Q X 3400 CONTINUE C. K=KM X I = IJL + (KM-1)*KK X SMJ = SMU*DJ(I) X DO 3600 N = 1, 4 X D4Q = Q(N,I-KK-KK) - 4.*(Q(N,I-KK)+Q(N,I+KK)) + 7.*Q(N,I) X IF (N .EQ. 3) D4Q = D4Q - 2.*Q(3,I) X S(N,I) = S(N,I) - SMJ*D4Q X 3600 CONTINUE X GO TO 5200 X 3800 CONTINUE X DO 5000 NK = 1, 2 X GO TO (4000,4200) NK X 4000 CONTINUE X KBGN = 1 X KEND = 2 X GO TO 4400 X 4200 CONTINUE X KBGN = KM X KEND = KMAX X 4400 CONTINUE X DO 4800 K = KBGN, KEND X KR = KKR(K) X KP = KKP(K) X KRR = KKR(KR) X KPP = KKP(KP) X I = IJL + (K-1)*KK X IR = IJL + (KR-1)*KK X IP = IJL + (KP-1)*KK X IRR = IJL + (KRR-1)*KK X IPP = IJL + (KPP-1)*KK X SMJ = SMU*DJ(I) X DO 4600 N = 1, 4 X D4Q=Q(N,IRR)+6.*Q(N,I)+Q(N,IPP)-4.*(Q(N,IR)+Q(N,IP)) X S(N,I) = S(N,I) - SMJ*D4Q X 4600 CONTINUE X 4800 CONTINUE X 5000 CONTINUE X 5200 CONTINUE X 5400 CONTINUE X 5600 CONTINUE X 5800 CONTINUE X DO 7000 J = 2, JM X DO 6800 K = KEND2, KENDM X IJK = J + (K-1)*KK X DO 6200 N = 1, 4 X L = 2 X I = IJK + (L-1)*LL X DO 6000 L = 3, LMM X I = I + LL X SMJ = SMU*DJ(I) X D4Q = Q(N,I-LL-LL) - 4.*(Q(N,I-LL)+Q(N,I+LL)) + 6.*Q(N,I) X 1 + Q(N,I+LL+LL) X S(N,I) = S(N,I) - SMJ*D4Q X 6000 CONTINUE X 6200 CONTINUE X I = IJK + LL X SMJ = SMU*DJ(I) X DO 6400 N = 1, 4 X D4Q = Q(N,I-LL) - 4.*Q(N,I) + 6.*Q(N,I+LL) - 4.*Q(N,I+LL+LL) X 1 + Q(N,I+LL+LL+LL) X S(N,I) = S(N,I) - SMJ*D4Q X 6400 CONTINUE X I = IJK + (LM-1)*LL X SMJ = SMU*DJ(I) X DO 6600 N = 1, 4 X D2Q = (-Q(N,I+LL)) + 2.*Q(N,I) - Q(N,I-LL) X S(N,I) = S(N,I) - SMJ*D2Q X 6600 CONTINUE X 6800 CONTINUE X 7000 CONTINUE X 7200 CONTINUE X RETURN X END X SUBROUTINE SUB021(NOP,NATOMS,NION,NSS) X COMMON /DATA/ FX(2000),FY(2000),FZ(2000) X COMMON /DATA/ X0(2000),Y0(2000),Z0(2000) X COMMON /DATA/ TX(2000),TY(2000),TZ(2000) X COMMON /DATA/ SX(2000),SY(2000),SZ(2000) X COMMON /DATA/ FSX(2000),FSY(2000),FSZ(2000) X COMMON /DATA/A(3),B(3),C(3),AA(3,650),BB(3,650),CC(3,650) X X ,AAA(650),BBB(650),CCC(650),AAAA,BBBB,CCCC X COMMON/DATA/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK X X ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3) X X ,R2(3),EWW,EWI,EWA,EII,EIA X COMMON/DATA/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30) X X ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP C COMMON/DATA/ QQ(11,11),AD(5),BD(5) X COMMON/DATA/ QQ(11,11), X 1A1,A2,A3,A4,A5,B1,B2,B3,B4,B5 X COMMON /DATA/ AA1(4),AA2(4),BB1(4),BB2(4),SFX(4),SFY(4),SFZ(4) X COMMON /DATA/ FOX(2000),FOY(2000),FOZ(2000), X 1F1X(2000),F1Y(2000), X 1F1Z(2000),F2X(2000),F2Y(2000),F2Z(2000), X 2FPX(2000),FPY(2000),FPZ(2000), X 3FAX(2000),FAY(2000),FAZ(2000),FIX(100),FIY(100),FIZ(100), X 4FOXP(2000),FOYP(2000),FOZP(2000),F1XP(2000),F1YP(2000), X 1F1ZP(2000),F2XP(2000),F2YP(2000),F2ZP(2000), X 2FPXP(2000),FPYP(2000),FPZP(2000), X 4XDT(2000),YDT(2000),ZDT(2000) X X DIMENSION NSITES(3),NSPECI(3),IND(2000) X DATA NSITES/3*1/,NSPECI/3*100/,IND/2000*0/ X X ALENGT = (22./UNITL)*.001 X ALENGM = 1./ALENGT X B1M = -B1 X B2M = -B2 X B3M = -B3 X B4M = -B4 X EFFE = 0 X EWW = 0 X EWI = 0 X EWA = 0 X EII = 0 X EIA = 0 X NSP = NSPECI(1) X NSST = NSS + NION + NATOMS X DO 1000 I = 1, NSST X FSX(I) = 0 X FSY(I) = 0 X FSZ(I) = 0 X 1000 CONTINUE X ISIT = NSITES(1) X RCUTS = 64.E-20/UNITL/UNITL X DO 1200 I = 1, NSP X FPX(I) = 0 X FPY(I) = 0 X FPZ(I) = 0 X FOX(I) = 0 X FOY(I) = 0 X FOZ(I) = 0 X F1X(I) = 0 X F1Y(I) = 0 X F1Z(I) = 0 X F2X(I) = 0 X F2Y(I) = 0 X F2Z(I) = 0 X 1200 CONTINUE X DO 2800 I = 2, NSP X INS = (I-1)*ISIT X FXO = 0 X FYO = 0 X FZO = 0 X FX1 = 0 X FY1 = 0 X FZ1 = 0 X FX2 = 0 X FY2 = 0 X FZ2 = 0 X FXP = 0 X FYP = 0 X FZP = 0 X DO 1600 J = 1, I - 1 X IND(J) = 0 X JNS = (J-1)*ISIT X XD = X0(I) - X0(J) X YD = Y0(I) - Y0(J) X ZD = Z0(I) - Z0(J) X XDT(J) = XD - 2.*ALENGT*(XD*ALENGM) X YDT(J) = YD - 2.*ALENGT*(YD*ALENGM) X ZDT(J) = ZD - 2.*ZD X DXS = XDT(J) + SX(INS+1) X DYS = YDT(J) + SY(INS+1) X DZS = ZDT(J) + SZ(INS+1) X RX = DXS - SX(JNS+1) X RY = DYS - SY(JNS+1) X RZ = DZS - SZ(JNS+1) X RSQ = RX*RX + RY*RY + RZ*RZ X IF (RSQ .GE. RCUTS) GO TO 1400 X IND(J) = 1 X 1400 CONTINUE X 1600 CONTINUE X L = 0 X DO 2000 J = 1, I - 1 X IF (IND(J) .EQ. 0) GO TO 1800 X L = L + 1 X IND(L) = J X 1800 CONTINUE X 2000 CONTINUE X IF (L .EQ. 0) GO TO 2600 X DO 2200 J = 1, L X K = IND(J) X JNS = (K-1)*ISIT X XD = XDT(K) X YD = YDT(K) X ZD = ZDT(K) X DXS = XD + SX(INS+1) X DYS = YD + SY(INS+1) X DZS = ZD + SZ(INS+1) X RX = DXS - SX(JNS+1) X RY = DYS - SY(JNS+1) X RZ = DZS - SZ(JNS+1) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X EWWT = A1*EXP(B1M*R) X FOR = (B1*EWWT)*R X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FXOT = FFX X FYOT = FFY X FZOT = FFZ X FOXJ = FFX X FOYJ = FFY X FOZJ = FFZ X RX = DXS - SX(JNS+2) X RY = DYS - SY(JNS+2) X RZ = DZS - SZ(JNS+2) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X ETERM1 = A2*EXP(B2M*R) X ETERM2 = A3*EXP(B3M*R) X FOR = (B2*ETERM1+B3*ETERM2)*R X EWWT = EWWT + ETERM1 + ETERM2 X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FXOT = FXOT + FFX X FYOT = FYOT + FFY X FZOT = FZOT + FFZ X F1XJ = FFX X F1YJ = FFY X F1ZJ = FFZ X RX = DXS - SX(JNS+3) X RY = DYS - SY(JNS+3) X RZ = DZS - SZ(JNS+3) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X ETERM1 = A2*EXP(B2M*R) X ETERM2 = A3*EXP(B3M*R) X FOR = (B2*ETERM1+B3*ETERM2)*R X EWWT = EWWT + ETERM1 + ETERM2 X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FXOT = FXOT + FFX X FYOT = FYOT + FFY X FZOT = FZOT + FFZ X F2XJ = FFX X F2YJ = FFY X F2ZJ = FFZ X DXS = XD + SX(INS+2) X DYS = YD + SY(INS+2) X DZS = ZD + SZ(INS+2) X RX = DXS - SX(JNS+2) X RY = DYS - SY(JNS+2) X RZ = DZS - SZ(JNS+2) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(2,3)*RM1 X EFFET = TEX X FOR = TEX*RSQ X ETERM1 = A4*EXP(B4M*R) X FOR = (B4*ETERM1)*RM1 + FOR X EWWT = EWWT + ETERM1 + TEX X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FX1T = FFX X FY1T = FFY X FZ1T = FFZ X F1XJ = F1XJ + FFX X F1YJ = F1YJ + FFY X F1ZJ = F1ZJ + FFZ X RX = DXS - SX(JNS+3) X RY = DYS - SY(JNS+3) X RZ = DZS - SZ(JNS+3) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(2,3)*RM1 X EFFET = EFFET + TEX X FOR = TEX*RSQ X ETERM1 = A4*EXP(B4M*R) X FOR = (B4*ETERM1)*RM1 + FOR X EWWT = EWWT + ETERM1 + TEX X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FX1T = FX1T + FFX X FY1T = FY1T + FFY X FZ1T = FZ1T + FFZ X F2XJ = F2XJ + FFX X F2YJ = F2YJ + FFY X F2ZJ = F2ZJ + FFZ X RX = DXS - SX(JNS+1) X RY = DYS - SY(JNS+1) X RZ = DZS - SZ(JNS+1) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X ETERM1 = A2*EXP(B2M*R) X ETERM2 = A3*EXP(B3M*R) X FOR = (B2*ETERM1+B3*ETERM2)*R X EWWT = EWWT + ETERM1 + ETERM2 X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FX1T = FX1T + FFX X FY1T = FY1T + FFY X FZ1T = FZ1T + FFZ X FOXJ = FOXJ + FFX X FOYJ = FOYJ + FFY X FOZJ = FOZJ + FFZ X RX = DXS - SX(JNS+4) X RY = DYS - SY(JNS+4) X RZ = DZS - SZ(JNS+4) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(2,4)*RM1 X EFFET = EFFET + TEX X FOR = TEX*RSQ X EWWT = EWWT + TEX X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FX1T = FX1T + FFX X FY1T = FY1T + FFY X FZ1T = FZ1T + FFZ X FPXJ = FFX X FPYJ = FFY X FPZJ = FFZ X DXS = XD + SX(INS+3) X DYS = YD + SY(INS+3) X DZS = ZD + SZ(INS+3) X RX = DXS - SX(JNS+2) X RY = DYS - SY(JNS+2) X RZ = DZS - SZ(JNS+2) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(2,3)*RM1 X EFFET = EFFET + TEX X FOR = TEX*RSQ X ETERM1 = A4*EXP(B4M*R) X FOR = (B4*ETERM1)*RM1 + FOR X EWWT = EWWT + ETERM1 + TEX X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FX2T = FFX X FY2T = FFY X FZ2T = FFZ X F1XJ = F1XJ + FFX X F1YJ = F1YJ + FFY X F1ZJ = F1ZJ + FFZ X RX = DXS - SX(JNS+3) X RY = DYS - SY(JNS+3) X RZ = DZS - SZ(JNS+3) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(2,3)*RM1 X EFFET = EFFET + TEX X FOR = TEX*RSQ X ETERM1 = A4*EXP(B4M*R) X FOR = (B4*ETERM1)*RM1 + FOR X EWWT = EWWT + ETERM1 + TEX X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FX2T = FX2T + FFX X FY2T = FY2T + FFY X FZ2T = FZ2T + FFZ X F2XJ = F2XJ + FFX X F2YJ = F2YJ + FFY X F2ZJ = F2ZJ + FFZ X RX = DXS - SX(JNS+1) X RY = DYS - SY(JNS+1) X RZ = DZS - SZ(JNS+1) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X ETERM1 = A2*EXP(B2M*R) X ETERM2 = A3*EXP(B3M*R) X FOR = (B2*ETERM1+B3*ETERM2)*R X EWWT = EWWT + ETERM1 + ETERM2 X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FX2T = FX2T + FFX X FY2T = FY2T + FFY X FZ2T = FZ2T + FFZ X FOXJ = FOXJ + FFX X FOYJ = FOYJ + FFY X FOZJ = FOZJ + FFZ X RX = DXS - SX(JNS+4) X RY = DYS - SY(JNS+4) X RZ = DZS - SZ(JNS+4) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(2,4)*RM1 X EFFET = EFFET + TEX X FOR = TEX*RSQ X EWWT = EWWT + TEX X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FX2T = FX2T + FFX X FY2T = FY2T + FFY X FZ2T = FZ2T + FFZ X FPXJ = FPXJ + FFX X FPYJ = FPYJ + FFY X FPZJ = FPZJ + FFZ X DXS = XD + SX(INS+4) X DYS = YD + SY(INS+4) X DZS = ZD + SZ(INS+4) X RX = DXS - SX(JNS+2) X RY = DYS - SY(JNS+2) X RZ = DZS - SZ(JNS+2) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(2,4)*RM1 X EFFET = EFFET + TEX X FOR = TEX*RSQ X EWWT = EWWT + TEX X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FXPT = FFX X FYPT = FFY X FZPT = FFZ X F1XJ = F1XJ + FFX X F1YJ = F1YJ + FFY X F1ZJ = F1ZJ + FFZ X RX = DXS - SX(JNS+3) X RY = DYS - SY(JNS+3) X RZ = DZS - SZ(JNS+3) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(2,4)*RM1 X EFFET = EFFET + TEX X FOR = TEX*RSQ X EWWT = EWWT + TEX X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FXPT = FXPT + FFX X FYPT = FYPT + FFY X FZPT = FZPT + FFZ X F2XJ = F2XJ + FFX X F2YJ = F2YJ + FFY X F2ZJ = F2ZJ + FFZ X RX = DXS - SX(JNS+4) X RY = DYS - SY(JNS+4) X RZ = DZS - SZ(JNS+4) X RSQ = RX*RX + RY*RY + RZ*RZ X R = SQRT(RSQ) X RM1 = R X TEX = QQ(4,4)*RM1 X EFFET = EFFET + TEX X FOR = TEX*RSQ X EWWT = EWWT + TEX X EWW = EWWT + EWW X EFFE = EFFET + EFFE X FFX = RX*FOR X FFY = RY*FOR X FFZ = RZ*FOR X FXPT = FXPT + FFX X FYPT = FYPT + FFY X FZPT = FZPT + FFZ X FPXJ = FPXJ + FFX X FPYJ = FPYJ + FFY X FPZJ = FPZJ + FFZ X FXO = FXOT + FXO X FYO = FYOT + FYO X FZO = FZOT + FZO X FX1 = FX1T + FX1 X FY1 = FY1T + FY1 X FZ1 = FZ1T + FZ1 X FX2 = FX2T + FX2 X FY2 = FY2T + FY2 X FZ2 = FZ2T + FZ2 X FXP = FXPT + FXP X FYP = FYPT + FYP X FZP = FZPT + FZP X FOXP(J) = FOXJ X FOYP(J) = FOYJ X FOZP(J) = FOZJ X F1XP(J) = F1XJ X F1YP(J) = F1YJ X F1ZP(J) = F1ZJ X F2XP(J) = F2XJ X F2YP(J) = F2YJ X F2ZP(J) = F2ZJ X FPXP(J) = FPXJ X FPYP(J) = FPYJ X FPZP(J) = FPZJ X 2200 CONTINUE X DO 2400 J = 1, L X K = IND(J) X FOX(K) = FOX(K) - FOXP(J) X FOY(K) = FOY(K) - FOYP(J) X FOZ(K) = FOZ(K) - FOZP(J) X F1X(K) = F1X(K) - F1XP(J) X F1Y(K) = F1Y(K) - F1YP(J) X F1Z(K) = F1Z(K) - F1ZP(J) X F2X(K) = F2X(K) - F2XP(J) X F2Y(K) = F2Y(K) - F2YP(J) X F2Z(K) = F2Z(K) - F2ZP(J) X FPX(K) = FPX(K) - FPXP(J) X FPY(K) = FPY(K) - FPYP(J) X FPZ(K) = FPZ(K) - FPZP(J) X 2400 CONTINUE X FOX(I) = FOX(I) + FXO X FOY(I) = FOY(I) + FYO X FOZ(I) = FOZ(I) + FZO X F1X(I) = F1X(I) + FX1 X F1Y(I) = F1Y(I) + FY1 X F1Z(I) = F1Z(I) + FZ1 X F2X(I) = F2X(I) + FX2 X F2Y(I) = F2Y(I) + FY2 X F2Z(I) = F2Z(I) + FZ2 X FPX(I) = FPX(I) + FXP X FPY(I) = FPY(I) + FYP X FPZ(I) = FPZ(I) + FZP X 2600 CONTINUE X 2800 CONTINUE X ISIT = NSITES(1) X NSP = NSPECI(1) X IF (NION .LE. 0) GO TO 3400 X ION = 0 X DO 3200 ION = 1, NION X K = NSP + ION X JSP = NSS + ION X XX = 0 X YY = 0 X ZZ = 0 X DO 3000 I = 1, NSP X INS = (I-1)*ISIT X XD = X0(I) - X0(K) X YD = Y0(I) - Y0(K) X ZD = Z0(I) - Z0(K) X XD = XD - 2.*ALENGT*(XD*ALENGM) X YD = YD - 2.*ALENGT*(YD*ALENGM) X ZD = ZD - 2.*ZD X DX = XD + SX(INS+1) X DY = YD + SY(INS+1) X DZ = ZD + SZ(INS+1) X DD = DX*DX + DY*DY + DZ*DZ X DDI = DD X DSQI = SQRT(DDI) X D6I = DDI*DDI*DDI X D12I = D6I*D6I X T1 = A(1)*DSQI X T2 = B(1)*D6I X T3 = C(1)*D12I X EWIT = T1 + T2 + T3 X GGG = T1 + 6.0*T2 + 12.0*T3 X GGG = GGG*DDI X XO = DX*GGG X YO = DY*GGG X ZO = DZ*GGG X DX = XD + SX(INS+2) X DY = YD + SY(INS+2) X DZ = ZD + SZ(INS+2) X DD = DX*DX + DY*DY + DZ*DZ X DDI = DD X DSQI = SQRT(DDI) X D6I = DDI*DDI*DDI X D12I = D6I*D6I X T1 = A(2)*DSQI X T2 = B(2)*D6I X T3 = C(2)*D12I X EWIT = EWIT + T1 + T2 + T3 X GGG = T1 + 6.0*T2 + 12.0*T3 X GGG = GGG*DDI X X1 = DX*GGG X Y1 = DY*GGG X Z1 = DZ*GGG X DX = XD + SX(INS+3) X DY = YD + SY(INS+3) X DZ = ZD + SZ(INS+3) X DD = DX*DX + DY*DY + DZ*DZ X DDI = DD X DSQI = SQRT(DDI) X D6I = DDI*DDI*DDI X D12I = D6I*D6I X T1 = A(3)*DSQI X T2 = B(3)*D6I X T3 = C(3)*D12I X EWIT = EWIT + T1 + T2 + T3 X EWI = EWI + EWIT X GGG = T1 + 6.0*T2 + 12.0*T3 X GGG = GGG*DDI X X2 = DX*GGG X Y2 = DY*GGG X Z2 = DZ*GGG X FOX(I) = FOX(I) + XO X FOY(I) = FOY(I) + YO X FOZ(I) = FOZ(I) + ZO X F1X(I) = F1X(I) + X1 X F1Y(I) = F1Y(I) + Y1 X F1Z(I) = F1Z(I) + Z1 X F2X(I) = F2X(I) + X2 X F2Y(I) = F2Y(I) + Y2 X F2Z(I) = F2Z(I) + Z2 X XX = XX + XO + X1 + X2 X YY = YY + YO + Y1 + Y2 X ZZ = ZZ + ZO + Z1 + Z2 X 3000 CONTINUE X FSX(JSP) = FSX(JSP) - XX X FSY(JSP) = FSY(JSP) - YY X FSZ(JSP) = FSZ(JSP) - ZZ X 3200 CONTINUE X 3400 CONTINUE X DO 3600 IA = 1, NATOMS X FAX(IA) = 0 X FAY(IA) = 0 X FAZ(IA) = 0 X 3600 CONTINUE X ISIT = NSITES(1) X NSP = NSPECI(1) X DO 4000 I = 1, NSP X INS = (I-1)*ISIT X DO 3800 IA = 1, NATOMS X K = IA + NOP X XD = X0(I) - X0(K) X YD = Y0(I) - Y0(K) X ZD = Z0(I) - Z0(K) X XD = XD - 2.*ALENGT*(XD*ALENGM) X YD = YD - 2.*ALENGT*(YD*ALENGM) X ZD = ZD - 2.*ZD X DX = XD + SX(INS+1) X DY = YD + SY(INS+1) X DZ = ZD + SZ(INS+1) X DD = DX*DX + DY*DY + DZ*DZ X DDI = DD X DSQI = SQRT(DDI) X D6I = DDI*DDI*DDI X D12I = D6I*D6I X T1 = AA(1,IA)*DSQI X T2 = BB(1,IA)*D6I X T3 = CC(1,IA)*D12I X EWAT = T1 + T2 + T3 X GGG = T1 + 6.0*T2 + 12.0*T3 X GGG = GGG*DDI X FOX(I) = FOX(I) + DX*GGG X FOY(I) = FOY(I) + DY*GGG X FOZ(I) = FOZ(I) + DZ*GGG X FAX(IA) = FAX(IA) + DX*GGG X FAY(IA) = FAY(IA) + DY*GGG X FAZ(IA) = FAZ(IA) + DZ*GGG X DX = XD + SX(INS+2) X DY = YD + SY(INS+2) X DZ = ZD + SZ(INS+2) X DD = DX*DX + DY*DY + DZ*DZ X DDI = DD X DSQI = SQRT(DDI) X D6I = DDI*DDI*DDI X D12I = D6I*D6I X T1 = AA(2,IA)*DSQI X T2 = BB(2,IA)*D6I X T3 = CC(2,IA)*D12I X EWAT = EWAT + T1 + T2 + T3 X GGG = T1 + 6.0*T2 + 12.0*T3 X GGG = GGG*DDI X F1X(I) = F1X(I) + DX*GGG X F1Y(I) = F1Y(I) + DY*GGG X F1Z(I) = F1Z(I) + DZ*GGG X FAX(IA) = FAX(IA) + DX*GGG X FAY(IA) = FAY(IA) + DY*GGG X FAZ(IA) = FAZ(IA) + DZ*GGG X DX = XD + SX(INS+3) X DY = YD + SY(INS+3) X DZ = ZD + SZ(INS+3) X DD = DX*DX + DY*DY + DZ*DZ X DDI = DD X DSQI = SQRT(DDI) X D6I = DDI*DDI*DDI X D12I = D6I*D6I X T1 = AA(3,IA)*DSQI X T2 = BB(3,IA)*D6I X T3 = CC(3,IA)*D12I X EWAT = EWAT + T1 + T2 + T3 X GGG = T1 + 6.0*T2 + 12.0*T3 X GGG = GGG*DDI X F2X(I) = F2X(I) + DX*GGG X F2Y(I) = F2Y(I) + DY*GGG X F2Z(I) = F2Z(I) + DZ*GGG X FAX(IA) = FAX(IA) + DX*GGG X FAY(IA) = FAY(IA) + DY*GGG X FAZ(IA) = FAZ(IA) + DZ*GGG X EWA = EWA + EWAT X 3800 CONTINUE X 4000 CONTINUE X DO 4200 I = 1, NSP X INS = (I-1)*ISIT X FSX(INS+1) = FSX(INS+1) + FOX(I) X FSY(INS+1) = FSY(INS+1) + FOY(I) X FSZ(INS+1) = FSZ(INS+1) + FOZ(I) X FSX(INS+2) = FSX(INS+2) + F1X(I) X FSY(INS+2) = FSY(INS+2) + F1Y(I) X FSZ(INS+2) = FSZ(INS+2) + F1Z(I) X FSX(INS+3) = FSX(INS+3) + F2X(I) X FSY(INS+3) = FSY(INS+3) + F2Y(I) X FSZ(INS+3) = FSZ(INS+3) + F2Z(I) X FSX(INS+4) = FSX(INS+4) + FPX(I) X FSY(INS+4) = FSY(INS+4) + FPY(I) X FSZ(INS+4) = FSZ(INS+4) + FPZ(I) X 4200 CONTINUE X IF (NION .LE. 0) GO TO 5800 X DO 4600 I = 1, NION X K = I + NSP X ISP = NSS + I X XX = 0 X YY = 0 X ZZ = 0 X XD = X0(K) X YD = Y0(K) X ZD = Z0(K) X DO 4400 J = 1, NATOMS X L = NOP + J X DX = XD - X0(L) X DY = YD - Y0(L) X DZ = ZD - Z0(L) X DX = DX - 2.*ALENGT*(DX*ALENGM) X DY = DY - 2.*ALENGT*(DY*ALENGM) X DZ = DZ - 2.*DZ X DD = DX*DX + DY*DY + DZ*DZ X DDI = DD X DSQI = SQRT(DDI) X D6I = DDI*DDI*DDI X D12I = D6I*D6I X T1 = AAA(J)*DSQI X T2 = BBB(J)*D6I X T3 = CCC(J)*D12I X EIA = EIA + T1 + T2 + T3 X GGG = T1 + 6.0*T2 + 12.0*T3 X GGG = GGG*DDI X XX = XX + GGG*DX X YY = YY + GGG*DY X ZZ = ZZ + GGG*DZ X FAX(J) = FAX(J) + GGG*DX X FAY(J) = FAY(J) + GGG*DY X FAZ(J) = FAZ(J) + GGG*DZ X 4400 CONTINUE X FSX(ISP) = FSX(ISP) + XX X FSY(ISP) = FSY(ISP) + YY X FSZ(ISP) = FSZ(ISP) + ZZ X 4600 CONTINUE X DO 4800 IA = 1, NATOMS X JSP = NSS + NION + IA X FSX(JSP) = FSX(JSP) - FAX(IA) X FSY(JSP) = FSY(JSP) - FAY(IA) X FSZ(JSP) = FSZ(JSP) - FAZ(IA) X 4800 CONTINUE X IF (NION .EQ. 1) GO TO 5800 X DO 5000 I = 1, NION X FIX(I) = 0 X FIY(I) = 0 X FIZ(I) = 0 X 5000 CONTINUE X DO 5400 I = 2, NION X XX = 0 X YY = 0 X ZZ = 0 X K = I + NSP X XD = X0(K) X YD = Y0(K) X ZD = Z0(K) X DO 5200 J = 1, I - 1 X L = NSP + J X JSP = NSS + J X DX = XD - X0(L) X DY = YD - Y0(L) X DZ = ZD - Z0(L) X DX = DX - 2.*ALENGT*(DX*ALENGM) X DY = DY - 2.*ALENGT*(DY*ALENGM) X DZ = DZ - 2.*DZ X DD = DX*DX + DY*DY + DZ*DZ X DDI = DD X DSQI = SQRT(DDI) X D6I = DDI*DDI*DDI X D12I = D6I*D6I X T1 = AAAA*DSQI X T2 = BBBB*D6I X T3 = CCCC*D12I X EII = EII + T1 + T2 + T3 X GGG = T1 + 6.0*T2 + 12.0*T3 X GGG = GGG*DDI X XX = XX + GGG*DX X YY = YY + GGG*DY X ZZ = ZZ + GGG*DZ X FSX(JSP) = FSX(JSP) - GGG*DX X FSY(JSP) = FSY(JSP) - GGG*DY X FSZ(JSP) = FSZ(JSP) - GGG*DZ X 5200 CONTINUE X FIX(I) = FIX(I) + XX X FIY(I) = FIY(I) + YY X FIZ(I) = FIZ(I) + ZZ X 5400 CONTINUE X DO 5600 I = 1, NION X ISP = NSS + I X FSX(ISP) = FSX(ISP) + FIX(I) X FSY(ISP) = FSY(ISP) + FIY(I) X FSZ(ISP) = FSZ(ISP) + FIZ(I) X 5600 CONTINUE X 5800 CONTINUE X DO 6000 I = 1, NOP X FX(I) = 0 X FY(I) = 0 X FZ(I) = 0 X 6000 CONTINUE X DO 6200 I = 1, NSPECI(1) X ISP = (I-1)*NSITES(1) X FX(I)=FX(I)+FSX(ISP+1)+FSX(ISP+2)+FSX(ISP+3)+FSX(ISP+4) X FY(I)=FY(I)+FSY(ISP+1)+FSY(ISP+2)+FSY(ISP+3)+FSY(ISP+4) X FZ(I)=FZ(I)+FSZ(ISP+1)+FSZ(ISP+2)+FSZ(ISP+3)+FSZ(ISP+4) X TX(I) = 0 X TY(I) = 0 X TZ(I) = 0 X TX(I) = FSZ(ISP+1)*SY(ISP+1) - FSY(ISP+1)*SZ(ISP+1) + FSZ(ISP+2) X 1 *SY(ISP+2) - FSY(ISP+2)*SZ(ISP+2) + FSZ(ISP+3)*SY(ISP+3) - X 2 FSY(ISP+3)*SZ(ISP+3) + FSZ(ISP+4)*SY(ISP+4) - FSY(ISP+4)*SZ( X 3 ISP+4) X TY(I) = FSX(ISP+1)*SZ(ISP+1) - FSZ(ISP+1)*SX(ISP+1) + FSX(ISP+2) X 1 *SZ(ISP+2) - FSZ(ISP+2)*SX(ISP+2) + FSX(ISP+3)*SZ(ISP+3) - X 2 FSZ(ISP+3)*SX(ISP+3) + FSX(ISP+4)*SZ(ISP+4) - FSZ(ISP+4)*SX( X 3 ISP+4) X TZ(I) = FSY(ISP+1)*SX(ISP+1) - FSX(ISP+1)*SY(ISP+1) + FSY(ISP+2) X 1 *SX(ISP+2) - FSX(ISP+2)*SY(ISP+2) + FSY(ISP+3)*SX(ISP+3) - X 2 FSX(ISP+3)*SY(ISP+3) + FSY(ISP+4)*SX(ISP+4) - FSX(ISP+4)*SY( X 3 ISP+4) X 6200 CONTINUE X N = NSS X DO 6400 I = NSPECI(1) + 1, 2000 X N = N + 1 X FX(I) = FSX(N) X FY(I) = FSY(N) X FZ(I) = FSZ(N) X 6400 CONTINUE X RETURN X END X SUBROUTINE SUB022(BUF,IW,NBES,K) X PARAMETER (NX=125, NY=125, NZ=10) X X COMMON /DATA/ BES(150),DZR X COMMON /DATA/ PP1(NX,NY,NZ),PM1(NX,NY,NZ) X COMMON /DATA/ CDPP1(NX,NY,NZ),CDPM1(NX,NY,NZ) X X DIMENSION BUF(NX,NY,2,1) X X DO 1200 J = 1, NY X DO 1000 I = 1, NX X PP1(I,J,K) = BUF(I,J,2,IW) X CDPP1(I,J,K) = BUF(I,J,1,IW) X PM1(I,J,K) = BUF(I,J,1,IW) X CDPM1(I,J,K) = BUF(I,J,2,IW) X 1000 CONTINUE X 1200 CONTINUE X DO 1600 J = 1, NY X DO 1400 I = 1, NX X BUF(I,J,1,IW) = BES(1)*PM1(I,J,K) + BES(2)*PP1(I,J,K) X BUF(I,J,2,IW) = BES(1)*CDPM1(I,J,K) + BES(2)*CDPP1(I,J,K) X 1400 CONTINUE X 1600 CONTINUE X DZR = DZR*2. X DO 2200 NB = 3, NBES X DO 2000 J = 1, NY X DO 1800 I = 1, NX X BUF(I,J,1,IW) = BUF(I,J,1,IW) + BES(NB)*PM1(I,J,K) X BUF(I,J,2,IW) = BUF(I,J,2,IW) + BES(NB)*CDPM1(I,J,K) X CC = PM1(I,J,K) X CCC = CDPM1(I,J,K) X PM1(I,J,K) = PP1(I,J,K) X CDPM1(I,J,K) = CDPP1(I,J,K) X PP1(I,J,K) = CC X CDPP1(I,J,K) = CCC X 1800 CONTINUE X 2000 CONTINUE X 2200 CONTINUE X RETURN X END X SUBROUTINE SUB023(VLCP) X LOGICAL ORTHO X PARAMETER (EPS=1.0E-6,EXCUT=30.0) X PARAMETER (IUMX = 1) X COMPLEX STRF X PARAMETER ( MXH=4414, MXV=401, MPX=21, NPK= 60, NAX=17, NTYPX= 3, X & LCORX=2, NPSX=3, NRX1=33, NRX2=32, NRX3=32, MAXTER=20, X & NRXX=NRX1*NRX2*NRX3,MXHC= MXH/1.415,MXVC= MXV/1.415) X PARAMETER ( NGM=3000,NTYP=NTYPX,NAT=NAX,NLC=2,NNL=3) X PARAMETER ( DQ=0.034, NQX=127 ) X PARAMETER ( NQXE=NQX*(NQX+1)/2) X COMMON/DATA/ TAU(3,NAX), X + PI, TPI, FPI, ALAT, TPIBA, TPIBA2, OMEGA, X + G(3,NRXX), GG(NRXX), X + ECUT(MAXTER), X + ZP(NPSX), CC(2,NPSX), ALPC(2,NPSX), X + APS(6,0:3,NPSX), ALPS(3,0:3,NPSX), X + XP(NPSX,NTYPX), X + ZV(NTYPX), X + TAB(NQX*(NQX+1)/2,0:LCORX,NTYPX), X + STRF(NRXX,NTYPX), VLOC(NRXX,NTYPX) X X DIMENSION LMAX(NPSX), LCORE(NPSX),ITYP(NAX),