c c bf 4/8/93: c This file contains the following routines: c pdedif, opdifsv, ckinpt, mvtou, driveb, mvtoun, mvunus c (note trsolv is in the update source twodsodium.src) c c SUBROUTINE PDEDIF(FUNC,BNBDY,NPDE,NZ,TEND,EPS,ZPTS,ITYPE,MBLFT, 1 MBRHT,NBUG,T,DT,INDEX,USOLN,WKCM,IWCM,NWKCM,NIWCM) SAVE C . C VERSION OF PDEDIF PACKAGE DATED 02/27/80 C SOLVES A PARABOLIC SYSTEM OF EQUATIONS IN ONE SPACE DIMENSION C AND TIME. TWO TYPES OF EQUATIONS ARE ALLOWED, FOR ITYPE=1 C . C DU/DT = D(A(U)*DU/DZ)/DZ + D(B(U)*U)/DZ + C(U,UZ) C . C AND FOR ITYPE=2 C . C DU/DT = A(U,UZ)*D2U/D2Z + B(U,UZ)*DU/DZ +C(U,UZ) C . C THE COEFFICIENTS A AND B ARE VECTOR RATHER THAN MATRIX, THAT IS C THERE IS NO COUPLING OF THE EQUATIONS THRU THE Z DERIVATIVES. C ACTUALLY THE COEFFICIENTS DEPEND ON Z AND T AS WELL AS U, C THAT IS A(M)=A(M,Z,T,U1,...UN) WHERE M IS THE EQUATION NUMBER C THE COEFFICIENTS A, B, AND C ARE COMPUTED BY THE USER SUPPLIED C ROUTINE FUNC. THE BOUNDARY CONDITIONS CAN BE EITHER C DIRICHLET OR MIXED, THE DIRICHLET MUST BE LINEAR IN U. THE C TYPE OF THE BOUNDARY CONDITION CAN BE DIFFERENT FOR EACH C COMPONENT OF THE SOLUTION VECTOR. C THE BOUNDARY CONDITIONS HAVE THE FORM C U = BL(T) AT THE LEFT OR LOWER BOUNDARY AND C U = BR(T) AT THE RIGHT OR UPPER BOUNDARY C THE FOLLOWING FORM IS ALSO ALLOWED C DU/DZ = BL(M,UL,T) OR DU/DZ = BR(M,UR,T) C A SECOND ORDER FINITE DIFFERENCE SCHEME IS USED IN SPACE C AND THE METHOD OF LINES IS USED FOR TIME, THE GEARB C INTEGRATOR OF HINDMARSH IS USED FOR TIME INTEGRATION. C C THE FUNC SUBROUTINE COMPUTES THE COEFFICIENTS A,B, AND C AND C RETURNS THE VALUES IN THE ARRAYS AN,BN, AND CN . C . C SUBROUTINE FUNC(NPDE,NZP2,ZN,T,UN,UZN,AN,BN,CN) C . C HERE AN,BN, AND CN ARE OUTPUT PARAMETERS AND THE REMAINING C ARGUMENTS ARE INPUT PARAMETERS. C THE DIMENSIONS OF THE ARGUMENT ARRAYS ARE C REAL ZN(NZP2),UN(NPDE,NZP2),UZN(NPDE,NZP2),AN(NPDE,NZP2), C 1 BN(NPDE,NZP2),CN(NPDE,NZP2) C WHERE NZP2=NZ+2. C NOTE THAT THE ORIGINAL MESH ARRAY ZPTS IS SHIFTED AND TWO C FICTITIOUS POINTS AT I=1 AND I=NZ+2 ARE ADDED TO FORM THE MESH C ARRAY ZN WHICH IS USED INTERNALLY IN THE PDEDIF PACKAGE. C BEWARE... THE ZPTS AND ZN ARRAYS ARE NOT THE SAME. C THE ARRAYS AN AND BN CAN DEPEND ON THE SOLUTION UN(M,I) C THE SOLUTION UN IS EXTRAPOLATED TO THE FICTITIOUS POINTS. C THE VALUES OF AN AND BN MUST BE COMPUTED WITHIN THIS ROUTINE AT C THESE FICTITIOUS POINTS. THEREFORE AN AND BN CAN NOT DEPEND ON C THE DERIVATIVE UZN OF UN. HOWEVER CN NEED NOT BE COMPUTED AT C THE FICTITIOUS POINTS AND THEREFORE CAN DEPEND ON UZN. NOTE C THAT UZN IS AN INPUT PARAMTETER OF FUNC. UZN IS NOT C DEFINED AT THE FICTITIOUS POINTS I=1 AND I=NZP2=NZ+2. C IF THE BOUNDARY TYPE IS MBLFT=1 OR MBLFT=3, THEN ALTHOUGH THE C VALUE OF AN AND BN AT THE FICTITIOUS POINTS MUST BE SET C THEY HAVE NO EFFECT ON THE SOLUTION. C . C THE SECOND USER SUPPLIED ROUTINE BNBDY IS USED TO DEFINE THE C BOUNDARY CONDITIONS. THIS ROUTINE HAS THE ARGUMENTS C . C SUBROUTINE BNBDY(NPDE,T,ZL,ZR,UL,UR,BL,BR) C . C WHERE THE ARGUMENT ARRAYS HAVE THE DIMENSIONS C . C REAL UL(NPDE),UR(NPDE),BL(NPDE),BR(NPDE) C . C . C HERE T IS THE TIME. THE LOWER OR LEFT BOUNDARY IS LOCATED AT C Z=ZL AND THE RIGHT OR UPPER BOUNDARY IS AT Z=ZR. ZL AND ZR C ARE INPUT PARAMETERS. THE BOUNDARY CONDITIONS HAVE THE FORM C . C U(M,ZL,T)=BL(M,T) IF MBLFT(M)=1 C OR C DU/DZ = BL(M,UL,T) IF MBLFT(M)=2 C OR C DU/DT(M,ZL,T) = BL(M,UL,T) IF MBLFT(M)=3 C . C THE UPPER OR RIGHT BOUNDARY IS SIMILAR. C . C U(M,ZR,T) = BR(M,T) C OR C DU/DZ(M,ZR,T) = BR(M,UR,T) C OR C DU/DT(M,ZR,T) = BR(M,UR,T) C . C . C . C THE INPUT/OUTPUT PARAMETERS ARE THE FOLLOWING C . C FUNC - THE USER SUPPLIED SUBROUTINE TO DEFINE THE RIGHT SIDE. C SEE ABOVE FOR MORE DESCRIPTION C BNBDY - THE USER SUPPLIED SUBROUTINE TO SUPPLY THE BOUNDARY C CONDITIONS. SEE ABOVE. C NPDE - THE NUMBER OF EQUATIONS OR SOLUTION VECTOR COMPONENTS C NZ - THE NUMBER OF MESH POINTS IN Z C TEND - THE STOPPING TIME FOR THE INTEGRATION C EPS - THE ERROR TOLERANCE FOR THE GEARB INTEGRATOR C ZPTS - THE ARRAY OF MESH POINTS ZPTS(NZ) C ITYPE - SELECTS THE FORMAT OF THE SYSTEM OF EQUATIONS C =1 FOR CONSERVATION FORM D(A*DU/DZ)/DZ C =2 FOR NON CONSERVATION FORM A*D2U/D2Z + B*DU/DZ + C C MBLFT - ARRAY GIVING TYPE OF BOUNDARY CONDITION FOR M-TH C SOLUTION COMPONENT AT THE LOWER BOUNDARY. C MBLFT(M) = 1 FOR U(M,ZL,T) = BL(M,T) C MBLFT(M) = 2 FOR DU/DZ(M,ZL,T) = BL(M,UL,T) C MBLFT(M) = 3 FOR DU/DT(M,ZL,T) = BL(M,UL,T) C MBRHT - ARRAY GIVING TYPE OF BOUNDARY CONDITION AT UPPER C OR RIGHT BOUNDARY. C NBUG - SELECTS DEBUGGING PRINTING C =1 TO PRINT MESSAGE WHEN GEARB REJECTS A TIME STEP BECAUSE C OF A LARGE ERROR ESTIMATE. C =2 TO PRINT CERTAIN GEARB PARAMETERS ON EACH TIME C STEP C =4 TO PRINT DIAGNOSTIC INFORMATION IN THE DIFFN ROUTINE C =8 TO PRINT DIAGNOSTIC INFOMMATION IN THE DIFFUN ROUTINE C T - ON INPUT THIS IS THE INITIAL TIME FOR THE INTEGRATION. C ON OUTPUT THIS IS THE TIME REACHED, WHICH IS TEND FOR C A SUCCESSFUL INTEGRATION. C DT - ON INPUT THIS GIVES THE INITIAL GUESS FOR THE TIME STEP C IF IN DOUBT, UNDERESTIMATE THIS VALUE. ON OUTPUT C THIS GIVES THE VALUE OF THE FINAL TIME STEP USED. C INDEX - DETERMINES HOW THE TIME INTEGRATION WILL BE STARTED C AND ALSO IS USED AS THE ERROR FLAG ON OUTPUT. C =1 TO START A NEW INTEGRATION C =0 TO CONTINUE AN INTEGRATION C =0 ON OUTPUT IF THE INTEGRATION WAS SUCCESSFUL. C USOLN - AN ARRAY HOLDING THE INITIAL VALUE OF THE SOLUTION ON C INPUT AND THE FINAL VALUE ON OUTPUT. THE DIMENSION IS C USOLN(NPDE,NZ) C WKCM - REAL WORK ARRAY. C NWKCM - DIMENSION OF WKCM. MUST EXCEED (24+6*NPDE)*NPDE*NZP2 C +5*NZP2 + 4*NPDE WHERE NZP2=NZ+2 C IWCM - INTEGER WORK ARRAY C NIWCM DIMENSION OF IWCM. MUST EXCEED 2*NPDE C . C . REAL ZPTS(NZ),USOLN(NPDE,NZ), WKCM(NWKCM) INTEGER IWCM(NIWCM) INTEGER MBLFT(NPDE),MBRHT(NPDE) COMMON /PDVAR/ NPDEC,NZC,NBUGC,ITYPEC COMMON /SYSPAR/ UROUND,MIN,MOUT EXTERNAL FUNC,BNBDY C . C CHECK THE INPUT DATA C . CALL CKINPT(ZPTS,NZ,NPDE,T,TEND,DT,MBLFT,MBRHT,NWKCM,NIWCM, 1 INDEX) IF(INDEX .LT. -31)RETURN C . C SET VARIABLES IN COMMON C . NPDEC=NPDE NZC=NZ ITYPEC=ITYPE NBUGC=NBUG C . C MOVE INPUT VARIABLES TO WORK ARRAY FOR COMMUNICATION WITH C OTHER ROUTINES C . NZP2=NZ+2 LML=1 LMR=LML+NPDE LZN=1 LUN=LZN+NZP2 LU=LUN+NPDE*NZP2 LAL=5*NZP2+7*NPDE*NZP2+1 LAR=LAL+NPDE LBL=LAR+NPDE LBR=LBL+NPDE LY=LBR+NPDE LYMAX=LY+13*NPDE*NZP2 LERROR=LYMAX+NPDE*NZP2 LSAV1=LERROR+NPDE*NZP2 LSAV2=LSAV1+NPDE*NZP2 LPW=LSAV2+NPDE*NZP2 N0PW=(6*NPDE-2)*NPDE*NZP2 LIPIV=LPW+N0PW DO 40 M=1,NPDE IWCM(M+LML-1)=MBLFT(M) 40 IWCM(M+LMR-1)=MBRHT(M) C . C ADD FICTITIOUS POINTS TO THE MESH C . DO 50 I=1,NZ WKCM(LZN+I)=ZPTS(I) DO 45 M=1,NPDE IM=LUN+M-1+I*NPDE 45 WKCM(IM)=USOLN(M,I) 50 CONTINUE WKCM(LZN)=ZPTS(1)-(ZPTS(2)-ZPTS(1)) WKCM(LZN+NZP2-1)=ZPTS(NZ)+(ZPTS(NZ)-ZPTS(NZ-1)) C . C COUNT THE NUMBER OF UNKNOWNS IN THE DIFFERENTIAL SYSTEM, THAT C IS SUBSTRACT THE DIRICHLET BOUNDARY VALUES AND THE FICTITIOUS C POINTS C . NK=0 DO 60 M=1,NPDE IF(MBLFT(M) .NE. 1)NK=NK+1 IF(MBRHT(M) .NE. 1)NK=NK+1 60 CONTINUE NEQ=NK + NPDE*(NZ-2) C . C MOVE INITIAL VALUES TO THE ARRAY U FOR INTEGRATION BY GEARB C NOTE THAT DIRICHLET BOUNDARY VALUES ARE DROPPED FROM U. C . CALL MVTOU(NPDE,NZ,NEQ,USOLN,MBLFT,MBRHT,WKCM(LU)) C . C CALL THE GEARB ODE SOLVER C . MF=22 ML=2*NPDE-1 MU=ML TBGN=T TLAST=TEND CALL DRIVEB(NEQ,TBGN,DT,WKCM(LU),TLAST,EPS,MF,INDEX,ML,MU, 1 N0PW,WKCM(LY),WKCM(LYMAX),WKCM(LERROR),WKCM(LSAV1), 2 WKCM(LSAV2),WKCM(LPW),WKCM(LIPIV),WKCM,IWCM,NWKCM,NIWCM,NBUG, 1 FUNC,BNBDY) C . C MOVE THE FINAL SOLUTION TO THE ARRAY USOLN C . CALL MVTOUN(TLAST,NPDE,NZP2,NEQ,WKCM(LU),MBLFT,MBRHT, 1 WKCM(LBL),WKCM(LBR),WKCM(LUN),WKCM(LZN),BNBDY) CALL MVUNUS(NPDE,NZP2,NZ,WKCM(LUN),USOLN) T=TLAST RETURN END SUBROUTINE OPDIFSV(UT,DTT) SAVE PARAMETER(KMX=96,KMX1=95,MMX=46,MMX1=45) include "blank.h" COMMON/SOLVTR/TP(MMX),DA(MMX),SHI(MMX),ATP(MMX),DAZ(MMX), 1SHIZ(MMX),ATPZ(MMX),ATPZZ(MMX),BNTR(MMX),CNTR(MMX),BRTR,QQOP(MMX) 2,XLOP(MMX),SHA(KMX),SHAZ(KMX),SHII(KMX),DAA(KMX) COMMON/TRAVLH/RK(50),ALP(20),RB(20),RKM(60) COMMON/WINDTRF/UN(KMX),VN(KMX),WN(KMX),UI(KMX),VI(KMX),WI(KMX) 1,WGAM(KMX),DWGAM(KMX),UIT,VUT DIMENSION USOLN(1,MMX),MBLFT(1),MBRHT(1), 1ZPTS(KMX) DIMENSION WZAM(MMX),DWZAM(MMX) COMMON/COLFACT/COLFAC REAL WKCN(20000) INTEGER IWCN(200) EXTERNAL FUNC,BNBDY COMMON/SYSPAR/UROUND,MIN,MOUT DATA MIN,MOUT/5,6/,UROUND/7.E-15/ C SOLVER PARAMETERS TM=0. EPS=1.E-2 DDT=1.E-2 DZIT=1.0 DZIU=0.5 ITYPE=2 NBUG=0 INDEX=1 NWKCN=20000 NIWCN=200 NPDE=1 MBLFT(1)=1 MBRHT(1)=2 ZMAS=16.*1.66E-24 DIP=50. SI=SIN(DIP/57.295) SII=SI*SI CI=COS(DIP/57.295) CII=CI*CI M1=KMX-MMX DO 1 K=1,MMX M=M1+K TP(K)=0.5*(TE(M)+TI(M)) TR=0.5*(TI(M)+TN(M)) A=ALOG10(TR) COLFAC=1.5 DA(K)=3.02E+17*TP(K)/(XNO(M)*SQRT(TP(K))*(1.08-0.139*A+4.5E-3*A* 1A)*COLFAC+19.9*XNN2(M)+19.5*XNO2(M)) DA(K)=DA(K)*SII SHI(K)=ZMAS*GZ(M)/(2.*BOLTZ*TP(K)) SHII(M)=1./SHI(K) DAA(M)=DA(K) SHA(K)=SHT(M) ATP(K)=ALOG(TP(K)) VINP=0.73E-9*SQRT(TN(M)/1000.)*XNO(M)*COLFAC+0.69E-9*XNN2(M)+ 10.67E-9*XNO2(M) BGAUS=0.4 XMASO=16. WOI=9.57946E+3*BGAUS/XMASO RHOI=VINP/WOI COI=1./(1+RHOI*RHOI) UIO=COI*(RHOI*RHOI*UN(M)+RHOI*(-VN(M)*SI-WN(M)*CI)) 1+COI*(UIT+RHOI*VUT) VIO=COI*(RHOI*RHOI*VN(M)+RHOI*UN(M)*SI+VN(M)*CII-WN(M)*SI*CI) 1+COI*(VUT*SI-RHOI*UIT*SI) WIO=COI*(RHOI*RHOI*WN(M)+RHOI*UN(M)*CI-VN(M)*SI*CI+ 1WN(M)*SII)+COI*(VUT*CI-RHOI*UIT*CI) WZAM(K)=WIO 1 CONTINUE DO 2 K=2,MMX1 DAZ(K)=(DA(K+1)-DA(K-1))/DZIT SHIZ(K)=(SHI(K+1)-SHI(K-1))/DZIT ATPZ(K)=(ATP(K+1)-ATP(K-1))/DZIT ATPZZ(K)=(ATP(K+1)+ATP(K-1)-2.*ATP(K))/(DZIT*DZIT) SHAZ(K)=(SHA(K+1)-SHA(K-1))/DZIT DWZAM(K)=(WZAM(K+1)-WZAM(K-1))/DZIT 2 CONTINUE DAZ(1)=(DA(2)-DA(1))/DZIU ATPZ(1)=(ATP(2)-ATP(1))/DZIU SHIZ(1)=(SHI(2)-SHI(1))/DZIU SHAZ(1)=(SHA(2)-SHA(1))/DZIU DWZAM(1)=(WZAM(2)-WZAM(1))/DZIU ATPZZ(1)=ATPZZ(2) DAZ(MMX)=(DA(MMX)-DA(MMX-1))/DZIU ATPZ(MMX)=(ATP(MMX)-ATP(MMX-1))/DZIU SHIZ(MMX)=(SHI(MMX)-SHI(MMX-1))/DZIU SHAZ(MMX)=(SHA(MMX)-SHA(MMX-1))/DZIU DWZAM(MMX)=(WZAM(MMX)-WZAM(MMX-1))/DZIU ATPZZ(MMX)=ATPZZ(MMX-1) DO 4 K=1,MMX M=M1+K BNTR(K)=DA(K)*(ATPZ(K)/SHA(K)+SHI(K)) 1+DAZ(K)/SHA(K)-WZAM(K) CNTR(K)=DAZ(K) 1/SHA(K)*(ATPZ(K)/SHA(K)+SHI(K))+DA(K)*(ATPZZ(K) 2/(SHA(K)*SHA(K))+SHIZ(K)/SHA(K)-SHAZ(K)*ATPZ(K)/(SHA(K)**3)) 3-DWZAM(K)/SHA(K) CALL RATECOE(TI(M),TE(M),TN(M)) XLOP(K)=RK(1)*XNO2(M)+RK(2)*XNN2(M)+RK(9)*XNH2(M)+RK(10)*XNH2O(M) 1+RK(12)*XNH(M)+RK(22)*XN2D(M) QQOP(K)=RK(8)*XINP(M)*XNO(M)+RK(11)*XIHP(M)*XNO(M) 4 CONTINUE BRTR=-(ATPZ(MMX)/SHA(MMX)+SHI(MMX))+WZAM(MMX)/DA(MMX) TEND=DTT DO 3 K=1,MMX M=M1+K ZPTS(K)=ZP(M) USOLN(1,K) =XIOP(M) 847 CONTINUE 3 CONTINUE CALL PDEDIF(FUNC,BNBDY,NPDE,MMX,TEND,EPS,ZPTS,ITYPE,MBLFT, 1MBRHT,NBUG,TM,DDT,INDEX,USOLN,WKCN,IWCN,NWKCN,NIWCN) DO 24 K=1,MMX M=M1+K XIOP(M)=AMAX1(USOLN(1,K),1.E-20) 24 CONTINUE RETURN END SUBROUTINE CKINPT(ZPTS,NZ,NPDE,TBGN,TEND,DT,MBLFT,MBRHT, 1 NWKCM,NIWCM,INDEX) SAVE C . C CHECK THE INPUT PARAMETERS FOR ERRORS C . REAL ZPTS(NZ) INTEGER MBLFT(NPDE),MBRHT(NPDE) COMMON /SYSPAR/ UROUND,MIN,MOUT DATA NZSAV,NPDSAV,ZPSAV/-1,-1,-1./ C . IDX=0 IF(INDEX .GE. -1 .AND. INDEX .LE. 3)GOTO 40 WRITE (MOUT,35) INDEX 35 FORMAT(5H ****,7H INDEX=,I6,16H IS OUT OF RANGE) IDX=IDX-32 GOTO 230 C . 40 IF(NPDE .GE. 1)GOTO 50 IDX=IDX-32 WRITE (MOUT,45) 45 FORMAT(//5H ****,12H NPDE .LT. 1) GOTO 230 C . 50 IF(NZ .GE. 3)GOTO 60 IDX=IDX-64 WRITE (MOUT,55) 55 FORMAT(//5H ****,10H NZ .LT. 3) GOTO 230 C . C . C CHECK THAT PARAMETERS WERE NOT CHANGED ON CONTINUATION OF C INTEGRATION C . 60 IF(INDEX .EQ. 1)GOTO 138 IF(NZ .NE. NZSAV)GOTO 134 IF(NPDE .NE. NPDSAV)GOTO 134 NZD2=NZ/2 + 1 ZP=ZPTS(1)+ZPTS(NZD2)+ZPTS(NZ) IF(ABS(ZP-ZPSAV) .LE. 100.*UROUND)GOTO 140 134 IDX=IDX-128 WRITE (MOUT,135) 135 FORMAT(5H ****,17H NZ NPDE OR ZPTS , 1 38HCHANGED ON CONTINUATION OF INTEGRATION / 2 5X,31HSET INDEX=1 FOR NEW INTEGRATION ) GOTO 230 138 NZSAV=NZ NPDSAV=NPDE NZD2=NZ/2 + 1 ZPSAV=ZPTS(1)+ZPTS(NZD2)+ZPTS(NZ) C . C CHECK FOR MONOTONE ZPTS C . 140 I1=0 NZP2=NZ+2 DO 145 K=2,NZ IF(ZPTS(K) .LE. ZPTS(K-1))I1=1 145 CONTINUE IF(I1 .EQ. 0)GOTO 160 WRITE (MOUT,150) 150 FORMAT(//5H ****,30H ZPTS(I) NOT INCREASING WITH I) IDX=IDX-256 C . 160 IF(TBGN .LT. TEND)GOTO 170 IDX=IDX-512 WRITE (MOUT,165) 165 FORMAT(//5H ****,15H TBGN .GE. TEND) C . 170 IF(DT .GT. 0.)GOTO 180 IDX=IDX-1024 WRITE (MOUT,175) 175 FORMAT(//5H ****,11H DT .LE. 0.) C . 180 I1=0 DO 195 M=1,NPDE IF(MBLFT(M) .LT. 1 .OR. MBLFT(M) .GT. 3)I1=1 IF(MBRHT(M) .LT. 1 .OR. MBRHT(M) .GT. 3)I1=1 195 CONTINUE IF(I1 .EQ. 0)GOTO 210 IDX=IDX-2048 WRITE (MOUT,205) 205 FORMAT(//5H ****,28H MBLFT OR MBRHT OUT OF RANGE) C . 210 NLIM=(24+6*NPDE)*NPDE*NZP2 + 5*NZP2 + 4*NPDE IF(NWKCM .GE. NLIM)GOTO 220 IDX=IDX-4096 WRITE (MOUT,215) NLIM 215 FORMAT(//5H ****,18H NWKCM MUST EXCEED,I6) C . 220 IF(NIWCM .GE. 2*NPDE)GOTO 230 WRITE (MOUT,225) 225 FORMAT(5H ****,23HNIWCM MUST BE GE 2*NPDE) IDX=IDX-8192 C . 230 CONTINUE IF(IDX .LT. 0)INDEX=IDX RETURN END SUBROUTINE MVTOU(NPDE,NZ,NEQ,USOLN,MBLFT,MBRHT,U) SAVE C . C MOVE SOLUTION FROM ARRAY USOLN TO ARRAY U. DROPS DIRICHLET C BOUNDARY VALUES C . REAL U(NEQ),USOLN(NPDE,NZ) INTEGER MBLFT(NPDE),MBRHT(NPDE) C . NZM1=NZ-1 C . KT=0 DO 20 M=1,NPDE IF(MBLFT(M) .EQ. 1)GOTO 20 KT=KT+1 U(KT)=USOLN(M,1) 20 CONTINUE NQLFT=KT DO 40 M=1,NPDE KT=NQLFT+M DO 35 I=2,NZM1 U(KT)=USOLN(M,I) 35 KT=KT+NPDE 40 CONTINUE KT=NQLFT+NPDE*(NZ-2) DO 50 M=1,NPDE IF(MBRHT(M) .EQ. 1)GOTO 50 KT=KT+1 U(KT)=USOLN(M,NZ) 50 CONTINUE RETURN END SUBROUTINE DRIVEB (N, T0, H0, Y0, TOUT, EPS, MF, INDEX, ML, MU, 1 N0PW,Y,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,WKCM,IWCM,NWKCM, 2 NIWCM,NBUG,FUNC,BNBDY) SAVE C THIS IS MODIFIED FROM THE MARCH 19, 1975 VERSION OF C GEARB, A PACKAGE FOR THE SOLUTION OF THE INITIAL VALUE C PROBLEM FOR SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS, C DY/DT = F(Y,T), Y = (Y(1),Y(2),...,Y(N)). C GEARB IS A VARIANT OF THE GEAR PACKAGE TO BE USED WHEN C THE JACOBIAN MATRIX DF/DY HAS BANDED OR NEARLY BANDED FORM. C SUBROUTINE DRIVEB IS A DRIVER ROUTINE FOR THE GEARB PACKAGE. C . C REFERENCES C 1. A. C. HINDMARSH, GEAR.. ORDINARY DIFFERENTIAL EQUATION C SYSTEM SOLVER, UCID-30001 REV. 3, LAWRENCE LIVERMORE C LABORATORY, P.O.BOX 808, LIVERMORE, CA 94550, DEC. 1974. C . C 2. A. C. HINDMARSH, GEARB.. SOLUTION OF ORDINARY C DIFFERENTIAL EQUATIONS HAVING BANDED JACOBIAN, C UCID-30059 REV. 1, L.L.L., MARCH 1975. C . C----------------------------------------------------------------------- C DRIVEB IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND C IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR, STIFFB. C . C THE INPUT PARAMETERS ARE.. C N = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS. C N CAN BE REDUCED, BUT NEVER INCREASED, DURING A PROBLEM. C T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE C (USED ONLY ON FIRST CALL). C H0 = THE NEXT STEP SIZE IN T (USED FOR INPUT ONLY ON THE C FIRST CALL). C Y0 = A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF C Y (USED FOR INPUT ONLY ON FIRST CALL). C TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. C INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT C AND THE PACKAGE WILL INTRRPOLATE TO T = TOUT. C EPS = THE RELATIVE ERROR BOUND (USED ONLY ON THE C FIRST CALL, UNLESS INDEX = -1). SINGLE STEP ERROR C ESTIMATES DIVIDED BY YMAX(I) WILL BE KEPT LESS THAN C EPS IN ROOT-MEAN-SQUARE NORM (I.E. EUCLIDEAN NORM C DIVIDED BY SQRT(N) ). THE VECTOR YMAX OF WEIGHTS C IS COMPUTED IN DRIVEB. INITIALLY YMAX(I) IS C ABS(Y(I)), WITH A DEFAULT VALUE OF 1 IF Y(I) = 0 C INITIALLY. THEREAFTER, YMAX(I) IS THE LARGEST VALUE C OF ABS(Y(I)) SEEN SO FAR, OR THE INITIAL YMAX(I) IF C THAT IS LARGER. TO ALTER EITHER OF THESE, CHANGE THE C APPROPRIATE STATEMENTS IN THE DO-LOOPS ENDING AT C STATEMENTS 10 AND 70 BELOW. C MF = THE METHOD FLAG (USED ONLY ON FIRST CALL, UNLESS C INDEX = -1). ALLOWED VALUES ARE 10, 11, 12, 13, C 20, 21, 22, 23. MF HAS TWO DECIMAL DIGITS, METH C AND MITER (MF = 10*METH + MITER). C METH IS THE BASIC METHOD INDICATOR.. C METH = 1 MEANS THE ADAMS METHODS. C METH = 2 MEANS THE BACKWARD DIFFERENTIATION C FORMULAS (BDF), OR STIFF METHODS OF GEAR. C MITER IS THE ITERATION METHOD INDICATOR.. C MITER = 0 MEANS FUNCTIONAL ITERATION (NO PARTIAL C DERIVATIVES NEEDED). C MITER = 1 MEANS CHORD METHOD WITH ANALYTIC JACOBIAN. C FOR THIS USER SUPPLIES SUBROUTINE C PDB (SEE DESCRIPTION BELOW). C MITER = 2 MEANS CHORD METHOD WITH JACOBIAN CALCULATED C INTERNALLY BY FINITE DIFFERENCES. C MITER = 3 MEANS CHORD METHOD WITH JACOBIAN REPLACED C BY A DIAGONAL APPROXIMATION BASED ON A C DIRECTIONAL DERIVATIVE. C INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL, C WITH THE FOLLOWING VALUES AND MEANINGS.. C 1 THIS IS THE FIRST CALL FOR THIS PROBLEM. C 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM, C AND INTEGRATION IS TO CONTINUE. C -1 THIS IS NOT THE FIRST CALL FOR THE PROBLEM, C AND THE USER HAS RESET N, EPS, AND/OR MF. C 2 SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT C EXACTLY (NO INTRRPOLATION IS DONE). C ASSUMES TOUT .GE. THE CURRENT T. C 3 SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING C PROGRAM AFTER ONE STEP. TOUT IS IGNORED. C SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0, C IT NEED NOT BE RESET FOR NORMAL CONTINUATION. C ML,MU = THE WIDTHS OF THE LOWER AND UPPER PARTS, RESPECTIVELY, C OF THE BAND IN THE JACOBIAN MATRIX, NOT COUNTING THE C MAIN DIAGONAL. THE FULL BANDWIDTH IS ML + MU + 1. C . C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN. C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. C A CHANGE OF PARAMETERS WITH INDEX = -1 CAN BE MADE AFTER C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN. C . C THE OUTPUT PARAMETERS ARE.. C H0 = THE STEP SIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT. C Y0 = THE COMPUTED VALUES OF Y AT T = TOUT. C TOUT = THE OUTPUT VALUE OF T. IF INTEGRATION WAS SUCCESSFUL, C AND THE INPUT VALUE OF INDEX WAS NOT 3, TOUT IS C UNCHANGED FROM ITS INPUT VALUE. OTHERWISE, TOUT C IS THE CURRENT VALUE OF T TO WHICH INTEGRATION C HAS BEEN COMPLETED. C INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS, C WITH THE FOLLOWING VALUES AND MEANINGS.. C 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND. C -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE C ERROR TEST EVEN AFTER REDUCING H BY A FACTOR OF C 1.E10 FROM ITS INITIAL VALUE. C -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS C HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY C A TEST ON EPS. TOO MUCH ACCURACY HAS BEEN REQUESTED. C -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE C CORRECTOR CONVERGENCE EVEN AFTER REDUCING H BY A C FACTOR OF 1.E10 FROM ITS INITIAL VALUE. C -4 IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT C PARAMETERS. SEE PRINTED MESSAGE. C -5 INDEX WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF C PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT C WAS NOT BEYOND T. INTRRPOLATION TO T = TOUT WAS C PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN, C SIMPLY CALL AGAIN WITH INDEX = -1 AND A NEW TOUT. C . C IN ADDITION TO DRIVEB, THE FOLLOWING ROUTINES ARE PROVIDED IN C THE PACKAGE.. C INTRRP(TOUT,Y,N0,Y0) INTRRPOLATES TO GET THE OUTPUT VALUES C AT T = TOUT, FROM THE DATA IN THE Y ARRAY. C STIFFB(Y,N0) IS THE CORE INTEGRATOR ROUTINE. IT PERFORMS A C SINGLE STEP AND ASSOCIATED ERROR CONTROL. C COSET(METH,NQ,EL,TQ,MAXDER) SETS COEFFICIENTS FOR USE IN C THE CORE INTEGRATOR. C PSETB(Y,N0,CON,MITER,IER) COMPUTES AND PROCESSES THE JACOBIAN C MATRIX J = DF/DY. C DECB(N0,N,ML,MU,B,IP,IER) AND SOLB(N0,N,ML,MU,B,C,X,IP) ARE C USED TO SOLVE BANDED LINEAR SYSTEMS. C NOTE.. PSETB, DECB, AND SOLB ARE CALLED ONLY IF MITER = 1 OR 2. C . C THE FOLLOWING ROUTINES ARE TO BE SUPPLIED BY THE USER.. C DIFFUN(N,T,Y,YDOT) COMPUTES THE FUNCTION YDOT = F(Y,T), THE C RIGHT-HAND SIDE OF THE O.D.E. C HERE Y AND YDOT ARE VECTORS OF LENGTH N. C PDB(N,T,Y,PD,N0,ML,MU) COMPUTES THE N BY N JACOBIAN MATRIX OF C PARTIAL DERIVATIVES, AND STORES IT IN PD C AS AN N0 BY ML+MU+1 ARRAY. PD(I,J-I+ML+1) C IS TO BE SET TO THE PARTIAL DERIVATIVE OF C YDOT(I) WITH RESPECT TO Y(J). CALLED ONLY C IF MITER = 1. OTHERWISE A DUMMY ROUTINE C CAN BE SUBSTITUTED. C . C THE DIMENSIONS IN THE FOLLOWING DECLARATIONS ARE SET FOR A C MAXIMUM OF 100 EQUATIONS AND FOR 2*ML+MU .LE. 30. IF THESE LIMITS C ARE TO BE EXCEEDED, THE DIMENSIONS SHOULD BE INCREASED ACCORDINGLY. C THE DIMENSION OF PW BELOW MUST BE AT LEAST N*(2*ML+MU+1) IF MITER = C 1 OR 2, BUT CAN BE REDUCED TO N IF MITER = 3, OR TO 1 IF MITER = 0. C THE DIMENSIONS OF YMAX, ERROR, SAVE1, SAVE2, IPIV, AND THE FIRST C DIMENSION OF Y SHOULD ALL BE AT LEAST N. THE COLUMN LENGTH OF C THE Y ARRAY AS USED ELSEWHERE IS N0, NOT 100. THE ROW LENGTH OF Y C CAN BE REDUCED FROM 13 TO 6 IF METH = 2. C THE IPIV ARRAY IS USED ONLY IF MITER IS 1 OR 2. C . C THE COMMON BLOCK GEAR9 CAN BE ACCESSED EXTERNALLY BY THE USER C IF DESIRED. IT CONTAINS THE STEP SIZE LAST USED (SUCCESSFULLY), C THE ORDER LAST USED (SUCCESSFULLY), THE NUMBER OF STEPS TAKEN C SO FAR, THE NUMBER OF F EVALUATIONS (DIFFUN CALLS) SO FAR, C AND THE NUMBER OF JACOBIAN EVALUATIONS SO FAR. C . C IN THE FOLLOWING DATA STATEMENT, SET.. C UROUND = THE UNIT ROUNDOFF OF THE MACHINE, I.E. THE SMALLEST C POSITIVE U SUCH THAT 1. + U .NE. 1. ON THE MACHINE. C LOUT = THE LOGICAL UNIT NUMBER FOR THE OUTPUT OF MESSAGES C DURING THE INTEGRATION. C----------------------------------------------------------------------- DIMENSION Y0(N) DIMENSION Y(N,13),WKCM(NWKCM),IWCM(NIWCM) COMMON /GEAR1/ T,H,HMIN,HMAX,EPSC,UROUND,NC,MFC,KFLAG,JSTART DIMENSION YMAX(N),ERROR(N),SAVE1(N),SAVE2(N),PW(N0PW),IPIV(N) COMMON /GEAR8/ EPSJ,MLC,MUC,MW,NM1,N0ML,N0W COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE C . COMMON /SYSPAR/ URUND,LDUMMY,LOUT EXTERNAL FUNC,BNBDY C . UROUND=URUND C . IF (INDEX .EQ. 0) GO TO 20 IF (INDEX .EQ. 2) GO TO 25 IF (INDEX .EQ. -1) GO TO 30 IF (INDEX .EQ. 3) GO TO 40 IF (INDEX .NE. 1) GO TO 430 IF (EPS .LE. 0.) GO TO 400 IF (N .LE. 0) GO TO 410 IF ((T0-TOUT)*H0 .GE. 0.) GO TO 420 C----------------------------------------------------------------------- C IF INITIAL VALUES OF YMAX OTHER THAN THOSE SET BELOW ARE DESIRED, C THEY SHOULD BE SET HERE. ALL YMAX(I) MUST BE POSITIVE. C IF VALUES FOR HMIN OR HMAX, THE BOUNDS ON ABS(H), OTHER THAN C THOSE BELOW ARE DESIRED, THEY SHOULD BE SET BELOW. C----------------------------------------------------------------------- DO 10 I = 1,N YMAX(I) = ABS(Y0(I)) IF (YMAX(I) .EQ. 0.) YMAX(I) = 1. 10 Y(I,1) = Y0(I) NC = N T = T0 H = H0 IF ((T+H) .EQ. T) WRITE(LOUT,15) 15 FORMAT(35H WARNING.. T + H = T ON NEXT STEP.) HMIN = ABS(H0) HMAX = ABS(T0-TOUT)*10. EPSC = EPS MFC = MF JSTART = 0 N0 = N EPSJ = SQRT(UROUND) MLC = ML MUC = MU MW = ML + MU + 1 NM1 = N0 - 1 N0ML = N0*ML N0W = N0*MW NHCUT = 0 GO TO 50 C . C TOUTP IS THE PREVIOUS VALUE OF TOUT FOR USE IN HMAX.----------------- 20 HMAX = ABS(TOUT-TOUTP)*10. GO TO 80 C . 25 HMAX = ABS(TOUT-TOUTP)*10. IF ((T-TOUT)*H .GE. 0.) GO TO 500 GO TO 85 C . 30 IF ((T-TOUT)*H .GE. 0.) GO TO 440 JSTART = -1 NC = N EPSC = EPS MFC = MF C . 40 IF ((T+H) .EQ. T) WRITE(LOUT,15) C . 50 CALL STIFFB (Y, N0,N0PW,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV, 1 WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) C . IF(NBUG/2-2*(NBUG/4) .NE. 1)GOTO 1055 WRITE (LOUT,1050) T,NSTEP,H,NFE,NJE,NQUSED 1050 FORMAT(2X,2HT=,E11.4,7H NSTEP=,I4,3H H=,E10.3,8H F EVAL=, 1 I4,10H JAC EVAL=,I4,7H ORDER=,I2) 1055 CONTINUE C . KGO = 1 - KFLAG GO TO (60, 100, 200, 300), KGO C KFLAG = 0, -1, -2, -3 C . 60 CONTINUE C----------------------------------------------------------------------- C NORMAL RETURN FROM INTEGRATOR. C . C THE WEIGHTS YMAX(I) ARE UPDATED. IF DIFFERENT VALUES ARE DESIRED, C THEY SHOULD BE SET HERE. A TEST IS MADE FOR EPS BEING TOO SMALL C FOR THE MACHINE PRECISION. C . C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY C STEP SHOULD BE INSERTED HERE. C . C IF INDEX = 3, Y0 IS SET TO THE CURRENT Y VALUES ON RETURN. C IF INDEX = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF C ERROR), AND THEN THE CURRENT Y VALUES ARE PUT IN Y0 ON RETURN. C FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR C UNLESS TOUT HAS BEEN REACHED. THEN INTRRPOLATED VALUES OF Y ARE C COMPUTED AND STORED IN Y0 ON RETURN. C IF INTRRPOLATION IS NOT DESIRED, THE CALL TO INTRRP SHOULD BE C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520. C----------------------------------------------------------------------- D = 0. DO 70 I = 1,N AYI = ABS(Y(I,1)) YMAX(I) = AMAX1(YMAX(I), AYI) 70 D = D + (AYI/YMAX(I))**2 D = D*(UROUND/EPS)**2 IF (D .GT. FLOAT(N)) GO TO 250 IF (INDEX .EQ. 3) GO TO 500 IF (INDEX .EQ. 2) GO TO 85 80 IF ((T-TOUT)*H .LT. 0.) GO TO 40 CALL INTRRP (TOUT, Y, N0, Y0) GO TO 520 85 IF (((T+H)-TOUT)*H .LE. 0.) GO TO 40 IF (ABS(T-TOUT) .LE. 100.*UROUND*HMAX) GO TO 500 IF ((T-TOUT)*H .GE. 0.) GO TO 500 H = (TOUT - T)*(1. - 4.*UROUND) JSTART = -1 GO TO 40 C----------------------------------------------------------------------- C ON AN ERROR RETURN FROM INTEGRATOR, AN IMMEDIATE RETURN OCCURS IF C KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE. C TO RECOVER, H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10 C TIMES BEFORE GIVING UP. C----------------------------------------------------------------------- 100 IF(NBUG - 2*(NBUG/2) .NE. 1)GOTO 110 WRITE (LOUT,105) T 105 FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/ 1 38H ERROR TEST FAILED WITH ABS(H) = HMIN/) 110 IF (NHCUT .EQ. 10) GO TO 150 NHCUT = NHCUT + 1 HMIN = .1*HMIN H = .1*H IF(NBUG-2*(NBUG/2) .NE. 1)GOTO 1115 WRITE (LOUT,115) H 115 FORMAT(24H H HAS BEEN REDUCED TO ,E16.8, 1 26H AND STEP WILL BE RETRIED//) 1115 JSTART = -1 GO TO 40 C . 150 WRITE (LOUT,155) 155 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//) GO TO 500 C . 200 WRITE (LOUT,205) T,H 205 FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,5H H =, 1 E16.8/52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//) GO TO 500 C . 250 WRITE (LOUT,255) T 255 FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/ 1 56H EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/) KFLAG = -2 GO TO 500 C . 300 WRITE (LOUT,305) T 305 FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/ 1 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/) GO TO 110 C . 400 WRITE (LOUT,405) 405 FORMAT(//28H ILLEGAL INPUT.. EPS .LE. 0.//) INDEX = -4 RETURN C . 410 WRITE (LOUT,415) 415 FORMAT(//25H ILLEGAL INPUT.. N .LE. 0//) INDEX = -4 RETURN C . 420 WRITE (LOUT,425) 425 FORMAT(//36H ILLEGAL INPUT.. (T0-TOUT)*H .GE. 0.//) INDEX = -4 RETURN C . 430 WRITE (LOUT,435) INDEX 435 FORMAT(//24H ILLEGAL INPUT.. INDEX =,I5//) INDEX = -4 RETURN C . 440 WRITE(LOUT,445) T,TOUT,H 445 FORMAT(//44H INDEX = -1 ON INPUT WITH (T-TOUT)*H .GE. 0./ 1 4H T =,E16.8,9H TOUT =,E16.8,6H H =,E16.8/ 1 44H INTRRPOLATION WAS DONE AS ON NORMAL RETURN./ 2 41H DESIRED PARAMETER CHANGES WERE NOT MADE.) CALL INTRRP (TOUT, Y, N0, Y0) INDEX = -5 RETURN C . 500 TOUT = T DO 510 I = 1,N 510 Y0(I) = Y(I,1) 520 INDEX = KFLAG TOUTP = TOUT H0 = HUSED IF (KFLAG .NE. 0) H0 = H RETURN C----------------------- END OF SUBROUTINE DRIVEB ---------------------- END SUBROUTINE MVTOUN(T,NPDE,NZP2,NEQ,U,MBLFT,MBRHT, 1 BL,BR,UN,ZN,BNBDY) SAVE C . C MOVE THE SOLUTION FROM THE U ARRAY TO THE UN ARRAY AND ADD C THE BOUNDARY VALUES USING THE BOUNDARY VALUES FROM BNBDY C . REAL U(NEQ),UN(NPDE,NZP2),BL(NPDE),BR(NPDE) 1 ,ZN(NZP2) INTEGER MBLFT(NPDE),MBRHT(NPDE) C . NZ=NZP2-2 NZP1=NZP2-1 KT=0 C . DO 20 M=1,NPDE IF(MBLFT(M) .EQ. 1)GOTO 20 KT=KT+1 UN(M,2)=U(KT) 20 CONTINUE NQLFT=KT DO 40 M=1,NPDE KT=NQLFT+M DO 35 I=3,NZ UN(M,I)=U(KT) 35 KT=KT+NPDE 40 CONTINUE KT=NQLFT+NPDE*(NZ-2) DO 50 M=1,NPDE IF(MBRHT(M) .EQ. 1)GOTO 50 KT=KT+1 UN(M,NZP2-1)=U(KT) 50 CONTINUE C . C SET DIRICHLET POINTS AND ALSO SET FICTITIOUS POINTS C . CALL BNBDY(NPDE,T,ZN(2),ZN(NZP2-1),UN(1,2),UN(1,NZP1),BL,BR) DO 125 M=1,NPDE C . C LOWER BOUNARY C . IF(MBLFT(M) .GT. 1)GOTO 110 C . C DIRICHLET BOUNDARY C . UN(M,2)=BL(M) UN(M,1)=UN(M,2)+((ZN(1)-ZN(2))/(ZN(3)-ZN(2)))*(UN(M,3)-UN(M,2)) GOTO 115 C . C NEUMANN BOUNDARY C . 110 UN(M,1)=UN(M,3)-(ZN(3)-ZN(1))*BL(M) C . C UPPER BOUNDARY C . 115 IF(MBRHT(M) .GT. 1)GOTO 120 C . C DIRICHLET BOUNDARY C . UN(M,NZP2-1)=BR(M) UN(M,NZP2)=UN(M,NZP2-1)+((ZN(NZP2)-ZN(NZP1))/(ZN(NZP1)-ZN(NZ)))* 1 (UN(M,NZP1)-UN(M,NZ)) GOTO 125 C . C NEUMANN BOUNDARY C . 120 UN(M,NZP2)=UN(M,NZ)+(ZN(NZP2)-ZN(NZ))*BR(M) 125 CONTINUE RETURN END SUBROUTINE MVUNUS(NPDE,NZP2,NZ,UN,USOLN) SAVE C . C MOVE SOLUTION IN UN ARRAY TO USOLN ARRAY C . REAL UN(NPDE,NZP2),USOLN(NPDE,NZ) C . DO 20 M=1,NPDE DO 10 I=1,NZ 10 USOLN(M,I)=UN(M,I+1) 20 CONTINUE RETURN END SUBROUTINE STIFFB (Y, N0,N0PW,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV, 1 WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) SAVE DIMENSION Y(N0,13), WKCM(NWKCM),IWCM(NIWCM) COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,N,MF,KFLAG,JSTART DIMENSION YMAX(N0),ERROR(N0),SAVE1(N0),SAVE2(N0),PW(N0PW), 1 IPIV(N0) COMMON /GEAR8/ EPSJ,ML,MU,MW,NM1,N0ML,N0W COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE EXTERNAL FUNC,BNBDY C----------------------------------------------------------------------- C STIFFB PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS. C STIFFB IS A VERSION FOR BANDED FORM OF THE JACOBIAN MATRIX. C COMMUNICATION WITH STIFFB IS DONE WITH THE FOLLOWING VARIABLES.. C . C Y AN N0 BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES C AND THEIR SCALED DERIVATIVES. LMAX IS 13 FOR THE ADAMS C METHODS AND 6 FOR THE GEAR METHODS. LMAX - 1 = MAXDER C IS THE MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET. C Y(I,J+1) CONTAINS THE J-TH DERIVATIVE OF Y(I), SCALED BY C H**J/FACTORIAL(J) (J = 0,1,...,NQ). C N0 A CONSTANT INTEGER .GE. N, USED FOR DIMENSIONING PURPOSES. C T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN. C H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. C HMIN, THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEP SIZE C HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY C TIME, BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE. C EPS THE RELATIVE ERROR BOUND. SEE DESCRIPTION IN DRIVER. C UROUND THE UNIT ROUNDOFF OF THE MACHINE. C N THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS. C MF THE METHOD FLAG. SEE DESCRIPTION IN DRIVER. C KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. C 0 THE STEP WAS SUCCESFUL. C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED C WITH ABS(H) = HMIN. C -2 THE REQUESTED ERROR IS SMALLER THAN CAN C BE HANDLED FOR THIS PROBLEM. C -3 CORRECTOR CONVERGENCE COULD NOT BE C ACHIEVED FOR ABS(H) = HMIN. C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND C THE Y ARRAY ARE AS OF THE BEGINNING OF THE LAST C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED. C JSTART AN INTEGER USED ON INPUT AND OUTPUT. C ON INPUT, IT HAS THE FOLLOWING VALUES AND MEANINGS.. C 0 PERFORM THE FIRST STEP. C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST. C .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF C H, EPS, N, AND/OR MF. C ON EXIT, JSTART IS NQ, THE CURRENT ORDER OF THE METHOD. C YMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL C ERRORS IN Y ARE COMPARED. C ERROR AN ARRAY OF N ELEMENTS. ERROR(I)/TQ(2) IS THE ESTIMATED C ONE-STEP ERROR IN Y(I). C SAVE1, TWO ARRAYS OF WORKING STORAGE, C SAVE2 EACH OF LENGTH N. C PW A BLOCK OF LOCATIONS USED FOR PARTIAL DERIVATIVES IF C MITER IS NOT 0. SEE DESCRIPTION IN DRIVER. C IPIV AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT C INFORMATION IF MITER = 1 OR 2. C ML,MU THE LOWER AND UPPER HALF BANDWIDTHS, RESPECTIVELY, OF C THE JACOBIAN. SEE DESCRIPTION IN DRIVER. C----------------------------------------------------------------------- DIMENSION EL(13),TQ(4) DATA EL(2)/1./, OLDL0/1./ KFLAG = 0 TOLD = T IF (JSTART .GT. 0) GO TO 200 IF (JSTART .NE. 0) GO TO 120 C----------------------------------------------------------------------- C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT IS C CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2 C FOR THE NEXT INCREASE. C----------------------------------------------------------------------- CALL DIFFUN (N, T, Y, SAVE1,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) DO 110 I = 1,N 110 Y(I,2) = H*SAVE1(I) METH = MF/10 MITER = MF - 10*METH NQ = 1 L = 2 IDOUB = 3 RMAX = 1.E4 RC = 0. CRATE = 1. EPSOLD = EPS HOLD = H MFOLD = MF NOLD = N NSTEP = 0 NSTEPJ = 0 NFE = 1 NJE = 0 IRET = 3 GO TO 130 C----------------------------------------------------------------------- C IF THE CALLER HAS CHANGED METH, COSET IS CALLED TO SET C THE COEFFICIENTS OF THE METHOD. IF THE CALLER HAS CHANGED C N, EPS, OR METH, THE CONSTANTS E, EDN, EUP, AND BND MUST BE RESET. C E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS C TO TEST FOR INCREASING THE ORDER, EDN FOR DECREASING THE ORDER. C BND IS USED TO TEST FOR CONVERGENCE OF THE CORRECTOR ITERATES. C IF THE CALLER HAS CHANGED H, Y MUST BE RESCALED. C IF H OR METH HAS BEEN CHANGED, IDOUB IS RESET TO L + 1 TO PREVENT C FURTHER CHANGES IN H FOR THAT MANY STEPS. C----------------------------------------------------------------------- 120 IF (MF .EQ. MFOLD) GO TO 150 MEO = METH MIO = MITER METH = MF/10 MITER = MF - 10*METH MFOLD = MF IF (MITER .NE. MIO) IWEVAL = MITER IF (METH .EQ. MEO) GO TO 150 IDOUB = L + 1 IRET = 1 130 CALL COSET (METH, NQ, EL, TQ, MAXDER) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDL0 OLDL0 = EL(1) 140 FN = FLOAT(N) EDN = FN*(TQ(1)*EPS)**2 E = FN*(TQ(2)*EPS)**2 EUP = FN*(TQ(3)*EPS)**2 BND = FN*(TQ(4)*EPS)**2 GO TO (160, 170, 200), IRET 150 IF ((EPS .EQ. EPSOLD) .AND. (N .EQ. NOLD)) GO TO 160 EPSOLD = EPS NOLD = N IRET = 1 GO TO 140 160 IF (H .EQ. HOLD) GO TO 200 RH = H/HOLD H = HOLD IREDO = 3 GO TO 175 170 RH = AMAX1(RH,HMIN/ABS(H)) 175 RH = AMIN1(RH,HMAX/ABS(H),RMAX) R1 = 1. DO 180 J = 2,L R1 = R1*RH DO 180 I = 1,N 180 Y(I,J) = Y(I,J)*R1 H = H*RH RC = RC*RH IDOUB = L + 1 IF (IREDO .EQ. 0) GO TO 690 C----------------------------------------------------------------------- C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY C MULTIPLYING THE Y ARRAY BY THE PASCAL TRIANGLE MATRIX. C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1). C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, OR THE CALLER HAS C CHANGED MITER, IWEVAL IS SET TO MITER TO FORCE THE PARTIALS TO BE C UPDATED, IF PARTIALS ARE USED. IN ANY CASE, THE PARTIALS C ARE UPDATED AT LEAST EVERY 20-TH STEP. C----------------------------------------------------------------------- 200 IF (ABS(RC-1.) .GT. 0.3) IWEVAL = MITER IF (NSTEP .GE. NSTEPJ+20) IWEVAL = MITER T = T + H DO 210 J1 = 1,NQ DO 210 J2 = J1,NQ J = (NQ + J1) - J2 DO 210 I = 1,N 210 Y(I,J) = Y(I,J) + Y(I,J+1) C----------------------------------------------------------------------- C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS C MADE ON THE R.M.S. NORM OF EACH CORRECTION, USING BND, WHICH C IS DEPENDENT ON EPS. THE SUM OF THE CORRECTIONS IS ACCUMULATED C IN THE VECTOR ERROR(I). THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR C LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1. C----------------------------------------------------------------------- 220 DO 230 I = 1,N 230 ERROR(I) = 0. M = 0 CALL DIFFUN (N, T, Y, SAVE2,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) NFE = NFE + 1 IF (IWEVAL .LE. 0) GO TO 290 C----------------------------------------------------------------------- C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED BEFORE C STARTING THE CORRECTOR ITERATION. IWEVAL IS SET TO 0 AS AN C INDICATOR THAT THIS HAS BEEN DONE. IF MITER = 1 OR 2, P IS C COMPUTED AND PROCESSED IN PSETB. IF MITER = 3, THE MATRIX USED C IS P = I - H*EL(1)*D, WHERE D IS A DIAGONAL MATRIX. C----------------------------------------------------------------------- IWEVAL = 0 RC = 1. NJE = NJE + 1 NSTEPJ = NSTEP GO TO (250, 240, 260), MITER 240 NFE = NFE + MW 250 CON = -H*EL(1) CALL PSETB (Y, N0, CON, MITER, IER,N0PW,YMAX,SAVE1,SAVE2,PW,IPIV, 1 WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) IF (IER .NE. 0) GO TO 420 GO TO 350 260 R = EL(1)*.1 DO 270 I = 1,N 270 PW(I) = Y(I,1) + R*(H*SAVE2(I) - Y(I,2)) CALL DIFFUN (N, T, PW, SAVE1,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) NFE = NFE + 1 HL0 = H*EL(1) DO 280 I = 1,N R0 = H*SAVE2(I) - Y(I,2) PW(I) = 1. D = .1*R0 - H*(SAVE1(I) - SAVE2(I)) SAVE1(I) = 0. IF (ABS(R0) .LT. UROUND*YMAX(I)) GO TO 280 IF (ABS(D) .EQ. 0.) GO TO 420 PW(I) = .1*R0/D SAVE1(I) = PW(I)*R0 280 CONTINUE GO TO 370 290 IF (MITER .NE. 0) GO TO (350, 350, 310), MITER C----------------------------------------------------------------------- C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM C THE RESULT OF THE LAST DIFFUN CALL. C----------------------------------------------------------------------- D = 0. DO 300 I = 1,N R = H*SAVE2(I) - Y(I,2) D = D + ( (R-ERROR(I))/YMAX(I) )**2 SAVE1(I) = Y(I,1) + EL(1)*R 300 ERROR(I) = R GO TO 400 C----------------------------------------------------------------------- C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR, C F SUB (M), AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND C SIDE AND P AS COEFFICIENT MATRIX, USING THE LU DECOMPOSITION C IF MITER = 1 OR 2. IF MITER = 3, THE COEFFICIENT H*EL(1) C IN P IS UPDATED. C----------------------------------------------------------------------- 310 PHL0 = HL0 HL0 = H*EL(1) IF (HL0 .EQ. PHL0) GO TO 330 R = HL0/PHL0 DO 320 I = 1,N D = 1. - R*(1. - 1./PW(I)) IF (ABS(D) .EQ. 0.) GO TO 440 320 PW(I) = 1./D 330 DO 340 I = 1,N 340 SAVE1(I) = PW(I)*(H*SAVE2(I) - (Y(I,2) + ERROR(I))) GO TO 370 350 DO 360 I = 1,N 360 SAVE1(I) = H*SAVE2(I) - (Y(I,2) + ERROR(I)) CALL SOLB (N0, N, ML, MU, PW, SAVE1, IPIV) 370 D = 0. DO 380 I = 1,N ERROR(I) = ERROR(I) + SAVE1(I) D = D + (SAVE1(I)/YMAX(I))**2 380 SAVE1(I) = Y(I,1) + EL(1)*ERROR(I) C----------------------------------------------------------------------- C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. C----------------------------------------------------------------------- 400 IF (M .NE. 0) CRATE = AMAX1(.9*CRATE,D/D1) IF ((D*AMIN1(1.,2.*CRATE)) .LE. BND) GO TO 450 D1 = D M = M + 1 IF (M .EQ. 3) GO TO 410 CALL DIFFUN (N, T, SAVE1, SAVE2,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) GO TO 290 C----------------------------------------------------------------------- C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. IF PARTIALS C ARE INVOLVED BUT ARE NOT UP TO DATE, THEY ARE REEVALUATED FOR THE C NEXT TRY. OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF NOT, A C NO-CONVERGENCE EXIT IS TAKEN. C----------------------------------------------------------------------- 410 NFE = NFE + 2 IF (IWEVAL .EQ. -1) GO TO 440 420 T = TOLD RMAX = 2. DO 430 J1 = 1,NQ DO 430 J2 = J1,NQ J = (NQ + J1) - J2 DO 430 I = 1,N 430 Y(I,J) = Y(I,J) - Y(I,J+1) IF (ABS(H) .LE. HMIN*1.00001) GO TO 680 RH = .25 IREDO = 1 GO TO 170 440 IWEVAL = MITER GO TO 220 C----------------------------------------------------------------------- C THE CORRECTOR HAS CONVERGED. IWEVAL IS SET TO -1 IF PARTIAL C DERIVATIVES WERE USED, TO SIGNAL THAT THEY MAY NEED UPDATING ON C SUBSEQUENT STEPS. THE ERROR TEST IS MADE AND CONTROL PASSES TO C STATEMENT 500 IF IT FAILS. C----------------------------------------------------------------------- 450 IF (MITER .NE. 0) IWEVAL = -1 NFE = NFE + M D = 0. DO 460 I = 1,N 460 D = D + (ERROR(I)/YMAX(I))**2 IF (D .GT. E) GO TO 500 C----------------------------------------------------------------------- C AFTER A SUCCESSFUL STEP, UPDATE THE Y ARRAY. C CONSIDER CHANGING H IF IDOUB = 1. OTHERWISE DECREASE IDOUB BY 1. C IF IDOUB IS THEN 1 AND NQ .LT. MAXDER, THEN ERROR IS SAVED FOR C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A C FACTOR OF AT LEAST 1.1. IF NOT, IDOUB IS SET TO 10 TO PREVENT C TESTING FOR THAT MANY STEPS. C----------------------------------------------------------------------- KFLAG = 0 IREDO = 0 NSTEP = NSTEP + 1 HUSED = H NQUSED = NQ DO 470 J = 1,L DO 470 I = 1,N 470 Y(I,J) = Y(I,J) + EL(J)*ERROR(I) IF (IDOUB .EQ. 1) GO TO 520 IDOUB = IDOUB - 1 IF (IDOUB .GT. 1) GO TO 700 IF (L .EQ. LMAX) GO TO 700 DO 490 I = 1,N 490 Y(I,LMAX) = ERROR(I) GO TO 700 C----------------------------------------------------------------------- C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES. C RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR C ONE LOWER ORDER. C----------------------------------------------------------------------- 500 KFLAG = KFLAG - 1 T = TOLD DO 510 J1 = 1,NQ DO 510 J2 = J1,NQ J = (NQ + J1) - J2 DO 510 I = 1,N 510 Y(I,J) = Y(I,J) - Y(I,J+1) RMAX = 2. IF (ABS(H) .LE. HMIN*1.00001) GO TO 660 IF (KFLAG .LE. -3) GO TO 640 IREDO = 2 PR3 = 1.E+20 GO TO 540 C----------------------------------------------------------------------- C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS C PR1, PR2, AND PR3 ARE COMPUTED, BY WHICH H COULD BE DIVIDED C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. C IN THE CASE OF FAILURE, PR3 = 1.E20 TO AVOID AN ORDER INCREASE. C THE SMALLEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE C ADDITIONAL SCALED DERIVATIVE. C----------------------------------------------------------------------- 520 PR3 = 1.E+20 IF (L .EQ. LMAX) GO TO 540 D1 = 0. DO 530 I = 1,N 530 D1 = D1 + ((ERROR(I) - Y(I,LMAX))/YMAX(I))**2 ENQ3 = .5/FLOAT(L+1) PR3 = ((D1/EUP)**ENQ3)*1.4 + 1.4E-6 540 ENQ2 = .5/FLOAT(L) PR2 = ((D/E)**ENQ2)*1.2 + 1.2E-6 PR1 = 1.E+20 IF (NQ .EQ. 1) GO TO 560 D = 0. DO 550 I = 1,N 550 D = D + (Y(I,L)/YMAX(I))**2 ENQ1 = .5/FLOAT(NQ) PR1 = ((D/EDN)**ENQ1)*1.3 + 1.3E-6 560 IF (PR2 .LE. PR3) GO TO 570 IF (PR3 .LT. PR1) GO TO 590 GO TO 580 570 IF (PR2 .GT. PR1) GO TO 580 NEWQ = NQ RH = 1./PR2 GO TO 620 580 NEWQ = NQ - 1 RH = 1./PR1 GO TO 620 590 NEWQ = L RH = 1./PR3 IF (RH .LT. 1.1) GO TO 610 DO 600 I = 1,N 600 Y(I,NEWQ+1) = ERROR(I)*EL(L)/FLOAT(L) GO TO 630 610 IDOUB = 10 GO TO 700 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1)) GO TO 610 C----------------------------------------------------------------------- C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. C IN ANY CASE H IS RESET ACCORDING TO RH AND THE Y ARRAY IS RESCALED. C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE. C----------------------------------------------------------------------- IF (NEWQ .EQ. NQ) GO TO 170 630 NQ = NEWQ L = NQ + 1 IRET = 2 GO TO 130 C----------------------------------------------------------------------- C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED. C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE C Y ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED. C AFTER A TOTAL OF 7 FAILURES, AN EXIT IS TAKEN WITH KFLAG = -2. C----------------------------------------------------------------------- 640 IF (KFLAG .EQ. -7) GO TO 670 RH = .1 RH = AMAX1(HMIN/ABS(H),RH) H = H*RH CALL DIFFUN (N, T, Y, SAVE1,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) NFE = NFE + 1 DO 650 I = 1,N 650 Y(I,2) = H*SAVE1(I) IWEVAL = MITER IDOUB = 10 IF (NQ .EQ. 1) GO TO 200 NQ = 1 L = 2 IRET = 3 GO TO 130 C----------------------------------------------------------------------- C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. C----------------------------------------------------------------------- 660 KFLAG = -1 GO TO 700 670 KFLAG = -2 GO TO 700 680 KFLAG = -3 GO TO 700 690 RMAX = 10. 700 HOLD = H JSTART = NQ RETURN C----------------------- END OF SUBROUTINE STIFFB ---------------------- END SUBROUTINE PSETB (Y, N0, CON, MITER, IER,N0PW,YMAX,SAVE1,SAVE2, 1 PW,IPIV,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) SAVE DIMENSION Y(N0,1),WKCM(NWKCM),IWCM(NIWCM) COMMON /GEAR1/ T,H,DUMMY(3),UROUND,N,IDUMMY(3) DIMENSION YMAX(N0),SAVE1(N0),SAVE2(N0),PW(N0PW),IPIV(N0) COMMON /GEAR8/ EPSJ,ML,MU,MW,NM1,N0ML,N0W EXTERNAL FUNC,BNBDY C----------------------------------------------------------------------- C PSETB IS CALLED BY STIFFB TO COMPUTE AND PROCESS THE MATRIX C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN. C THIS IS THE VERSION FOR BANDED FORM OF J. C J IS COMPUTED, EITHER BY THE USER-SUPPLIED ROUTINE PDB C IF MITER = 1, OR BY FINITE DIFFERENCING IF MITER = 2. C J IS STORED IN PW AND REPLACED BY P, USING CON = -H*EL(1). C THEN P IS SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR C LATER SOLUTION OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. C . C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION C WITH PSETB USES THE FOLLOWING.. C EPSJ = SQRT(UROUND), USED IN THE NUMERICAL JACOBIAN INCREMENTS. C MW = ML + MU + 1. C NM1 = N0 - 1. C N0ML = N0*ML. C N0W = N0*MW. C----------------------------------------------------------------------- IF (MITER .EQ. 2) GO TO 20 C IF MITER = 1, CALL PDB AND MULTIPLY BY SCALAR. ---------------------- CALL PDB (N, T, Y, PW, N0, ML, MU) DO 10 I = 1,N0W 10 PW(I) = PW(I)*CON GO TO 90 C IF MITER = 2, MAKE MW CALLS TO DIFFUN TO APPROXIMATE J. ------------- 20 D = 0. DO 30 I = 1,N 30 D = D + SAVE2(I)**2 R0 = ABS(H)*SQRT(D)*1.E03*UROUND DO 80 J = 1,MW KMAX = (N - J)/MW + 1 YM = 0. DO 40 K = 1,KMAX I = J + (K - 1)*MW 40 YM = AMAX1(YM,YMAX(I)) R = AMAX1(EPSJ*YM,R0) DO 50 K = 1,KMAX I = J + (K - 1)*MW 50 Y(I,1) = Y(I,1) + R CALL DIFFUN (N, T, Y, SAVE1,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY) D = CON/R J1 = J*N0 + N0ML DO 70 K = 1,KMAX JJ = J + (K - 1)*MW I1 = MAX0(JJ-MU,1) I2 = MIN0(JJ+ML,N) II = J1 - NM1*I1 DO 60 I = I1,I2 PW(II) = (SAVE1(I) - SAVE2(I))*D 60 II = II - NM1 Y(JJ,1) = Y(JJ,1) - R 70 J1 = J1 + N0W 80 CONTINUE C ADD IDENTITY MATRIX. ------------------------------------------------ 90 DO 100 I = 1,N 100 PW(N0ML+I) = PW(N0ML+I) + 1. C DO LU DECOMPOSITION ON P. ------------------------------------------- CALL DECB (N0, N, ML, MU, PW, IPIV, IER) RETURN C----------------------- END OF SUBROUTINE PSETB ----------------------- END SUBROUTINE PDB SAVE C . C DUMMY ROUTINE FOR COMPUTATION OF JACOBIAN. SHOULD NOT BE USED C . WRITE (MOUT,10) 10 FORMAT(5H ****,11HSTOP IN PDB) STOP END SUBROUTINE INTRRP (TOUT, Y, N0, Y0) SAVE DIMENSION Y0(N0),Y(N0,1) COMMON /GEAR1/ T,H,DUMMY(4),N,IDUMMY(2),JSTART C----------------------------------------------------------------------- C SUBROUTINE INTRRP COMPUTES INTRRPOLATED VALUES OF THE DEPENDENT C VARIABLE Y AND STORES THEM IN Y0. THE INTRRPOLATION IS TO THE C POINT T = TOUT, AND USES THE NORDSIECK HISTORY ARRAY Y, AS FOLLOWS.. C NQ C Y0(I) = SUM Y(I,J+1)*S**J , C J=0 C WHERE S = -(T-TOUT)/H. C----------------------------------------------------------------------- DO 10 I = 1,N 10 Y0(I) = Y(I,1) L = JSTART + 1 S = (TOUT - T)/H S1 = 1. DO 30 J = 2,L S1 = S1*S DO 20 I = 1,N 20 Y0(I) = Y0(I) + S1*Y(I,J) 30 CONTINUE RETURN C----------------------- END OF SUBROUTINE INTRRP ---------------------- END SUBROUTINE NADIFSV(UT,DTT) SAVE PARAMETER(KMX=96,MMX=96) include "blank.h" include "sodium.h" COMMON/TRAVLH/RK(50),ALP(20),RB(20),RKM(60) COMMON/COLFACT/COLFAC COMMON/TIMETRF/TIMEE c COMMON/SOLVT/TP(MMX),DA(MMX),SHI(MMX),DAZ(MMX), 1SHIZ(MMX),BNTR(MMX),CNTR(MMX),BRTR,QQOP(MMX) 2,XLOP(MMX),SHA(MMX),SHAZ(MMX),RR(MMX),SS(MMX),DAM(MMX), 3DAMZ(MMX),DEDY(MMX),DRZ(MMX),DSZ(MMX),ANTR(MMX) COMMON/EDDTRF/EDY(KMX) COMMON/WINDTRF/UN(KMX),VN(KMX),WN(KMX),UI(KMX),VI(KMX),WI(KMX) 1,WGAM(KMX),DWGAM(KMX),UIT,VUT DIMENSION USOLN(1,MMX),MBLFT(1),MBRHT(1), 1ZPTS(MMX) DIMENSION DTETI(KMX),DLNNE(KMX),DTN(KMX) REAL WKCN(20000) INTEGER IWCN(200) EXTERNAL NAFUNC,NABNBDY COMMON/SYSPAR/UROUND,MIN,MOUT DATA MIN,MOUT/5,6/,UROUND/7.E-15/ C SOLVER PARAMETERS XJNAP=1.7E-5*0.5 RKSI1=6.0E-10 RKSI2=8.0E-11 RKSI3=6.2E-10 RKC1=3.2E-31 RKCM1=1.4E-15 RKC2=3.2E-29 RKCM2=1.2E-17 RKC3=1.0E-9 RKC4=1.0E-9 RKC5=1.2E-6 RKC6=1.2E-6 RKSI4=0. RKSI5=0. RKS1=3.4E-10 RKS2=1.2E-11 RKS2S=2.8E-11 RKS2T=RKS2+RKS2S RKS3=1.9E-30 XJNAO2=1.E-4 XJNAOH=2.0E-3 RKS6=1.0E-13 RKS7=2.0E-10 RKS8=3.5E-11 RKS8S=1.0E-10 RKS8T=RKS8+RKS8S RKS9=1.0E-14 RKS10=1.0E-10 RKS11=1.0E-13 RKS12=4.0E-12 TM=0. EPS=1.E-4 DDT=1.E-10 DZIT=1.0 DZIU=0.5 ITYPE=2 NBUG=0 INDEX=1 NWKCN=20000 NIWCN=200 NPDE=1 MBLFT(1)=1 MBRHT(1)=2 ZMAS=23.*1.66E-24 M1=KMX-MMX DIP=50. SI=SIN(DIP/57.295) SII=SI*SI CI=COS(DIP/57.295) CII=CI*CI KMXM1=KMX-1 DO 11 K=2,KMXM1 DTETI(K)=((TE(K+1)+TI(K+1))-(TE(K-1)+TI(K-1)))/DZIT DLNNE(K)=(ALOG(XNE(K+1))-ALOG(XNE(K-1)))/DZIT DTN(K)=(TN(K+1)-TN(K-1))/DZIT 11 CONTINUE DTETI(1)=((TE(2)+TI(2))-(TE(1)+TI(1)))/DZIU DLNNE(1)=(ALOG(XNE(2))-ALOG(XNE(1)))/DZIU DTN(1)=(TN(2)-TN(1))/DZIU DTETI(KMX)=((TE(KMX)+TI(KMX))-(TE(KMXM1)+TI(KMXM1)))/DZIU DLNNE(KMX)=(ALOG(XNE(KMX))-ALOG(XNE(KMXM1)))/DZIU DTN(KMX)=(TN(KMX)-TN(KMXM1))/DZIU DO 1 K=1,MMX M=M1+K COLFAC=1.5 SHI(K)=ZMAS*GZ(M)/(BOLTZ*TI(K)) SHA(K)=SHT(M) RR(K)=SHI(K)+DTETI(K)/(SHA(K)*TI(K))+TE(K)*DLNNE(K)/(SHA(K)* 1TI(K)) SS(K)=DTN(K)/(SHA(K)*TN(K))+1./SHA(K) C VINP=0.73E-9*SQRT(TN(K)/1000.)*XNO(K)*COLFAC+0.69E-9*XNN2(K)+ C 10.67E-9*XNO2(K) VINP=3.E-10*(XNO(M)+XNO2(M)+XNN2(M)) XMASO=23. DAM(K)=BOLTZ*TI(K)/(XMASO*1.66E-24*VINP) DA(K)=DAM(K)*SII+EDY(K)*SHA(K)*SHA(K) BGAUS=0.4 WOI=9.57946E+3*BGAUS/XMASO RHOI=VINP/WOI COI=1./(1+RHOI*RHOI) UI(K)=COI*(RHOI*RHOI*UN(K)+RHOI*(-VN(K)*SI-WN(K)*CI)) 1+COI*(UIT+RHOI*VUT) VI(K)=COI*(RHOI*RHOI*VN(K)+RHOI*UN(K)*SI+VN(K)*CII-WN(K)*SI*CI) 1+COI*(VUT*SI-RHOI*UIT*SI) WI(K)=COI*(RHOI*RHOI*WN(K)+RHOI*UN(K)*CI-VN(K)*SI*CI+ 1WN(K)*SII)+COI*(VUT*CI-RHOI*UIT*CI) WGAM(K)=WI(K) 1 CONTINUE MMX1=MMX-1 DO 2 K=2,MMX1 DAZ(K)=(DA(K+1)-DA(K-1))/DZIT DAMZ(K)=(DAM(K+1)-DAM(K-1))/DZIT DEDY(K)=(EDY(K+1)-EDY(K-1))/DZIT DRZ(K)=(RR(K+1)-RR(K-1))/DZIT DSZ(K)=(SS(K+1)-RR(K-1))/DZIT SHAZ(K)=(SHA(K+1)-SHA(K-1))/DZIT DWGAM(K)=(WGAM(K+1)-WGAM(K-1))/DZIT 2 CONTINUE DAZ(1)=(DA(2)-DA(1))/DZIU DAMZ(1)=(DAM(2)-DAM(1))/DZIU DEDY(1)=(EDY(2)-EDY(1))/DZIU DRZ(1)=(RR(2)-RR(1))/DZIU DSZ(1)=(SS(2)-SS(1))/DZIU SHAZ(1)=(SHA(2)-SHA(1))/DZIU DWGAM(1)=(WGAM(2)-WGAM(1))/DZIU DAZ(MMX)=(DA(MMX)-DA(MMX-1))/DZIU DAMZ(MMX)=(DAM(MMX)-DAM(MMX-1))/DZIU DEDY(MMX)=(EDY(MMX)-EDY(MMX-1))/DZIU DRZ(MMX)=(RR(MMX)-RR(MMX-1))/DZIU DSZ(MMX)=(SS(MMX)-SS(MMX-1))/DZIU SHAZ(MMX)=(SHA(MMX)-SHA(MMX-1))/DZIU DWGAM(MMX)=(WGAM(MMX)-WGAM(MMX-1))/DZIU DO 4 K=1,MMX M=M1+K ANTR(K)=DA(K)/(SHA(K)*SHA(K)) BNTR(K)=DAZ(K)/(SHA(K)*SHA(K))-DA(K)*SHAZ(K)/(SHA(K)**3) 1+(DAM(K)*RR(K)+EDY(K)*SS(K)-WGAM(K))/SHA(K) CNTR(K)=(RR(K)*DAMZ(K)+SS(K)*DEDY(K)+DAM(K)*DRZ(K)+EDY(K)* 1DSZ(K)-DWGAM(K))/SHA(K) AS=0. ZPS=100. H1S=24. H2S=40. PNAP=AS*EXP(-((ZPHT(K)-ZPS)/H1S)**2)*EXP((ZPHT(K)-ZPS)/H2S) PNASP(K)=(XJNAP+RKSI1*XIO2P(K)+RKSI2*XINOP(K)+RKSI3*XIN2P(K)) 1*XNAS(K)+PNAP C IF(TIMEE.LT.6..OR.TIMEE.GT.18.) PNASP(K)=PNASP(K)*1.E-3 XM=XNO(K)+XNO2(K)+XNN2(K) ASP=RKC1*XNN2(K)*XM/(RKC3*XNCO2(K)+RKCM1*XM) BSP=RKC2*XNCO2(K)*XM CSP=RKCM2*XM+RKC4*XNH2O(K)+RKC5*XNE(K) XLNASP(K)=RKC1*XNN2(K)*XM+RKC2*XNCO2(K)*XM+RKSI4*XNO3(K) 1-RKCM1*XM*ASP-RKCM2*XM*(BSP+RKC3*XNCO2(K)*ASP)/CSP 4 CONTINUE BRTR=-(DAM(MMX)*RR(MMX)+EDY(MMX)*SS(MMX)-WGAM(MMX))*SHA(MMX)/ 1(DAM(MMX)+EDY(MMX)) TEND=DTT DO 3 K=1,MMX M=M1+K ZPTS(K)=ZP(M) USOLN(1,K) =XNASP(M) 847 CONTINUE 3 CONTINUE CALL PDEDIF(NAFUNC,NABNBDY,NPDE,MMX,TEND,EPS,ZPTS,ITYPE, 1MBLFT,MBRHT,NBUG,TM,DDT,INDEX,USOLN,WKCN,IWCN,NWKCN,NIWCN) DO 24 K=1,MMX M=M1+K C XNASP(M)=USOLN(1,K) IF (ZPHT(M).GT.100..AND.XNASP(M).LT.1) XNASP(M)=1. C XNAOP(K)=RKSI4*XNO3(K)*XNASP(K)/(RKSI5*XNE(K)) XNAPN2(K)=RKC1*XNN2(K)*XM*XNASP(K)/(RKC3*XNCO2(K)+RKCM1*XM) XNAPCO2(K)=(RKC2*XNCO2(K)*XM*XNASP(K)+RKC3*XNCO2(K)*XNAPN2(K)) 1/(RKCM2*XM+RKC4*XNH2O(K)+RKC5*XNE(K)) XNAPH2O(K)=RKC4*XNH2O(K)*XNAPCO2(K)/(RKC6*XNE(K)) 24 CONTINUE RETURN END SUBROUTINE NABNBDY(NPDE,TM,ZL,ZR,AL,AR,BL,BR) SAVE PARAMETER(KMX=96,MMX=96) include "blank.h" DIMENSION AL(NPDE),AR(NPDE),BL(NPDE),BR(NPDE) COMMON/SOLVT/TP(MMX),DA(MMX),SHI(MMX),DAZ(MMX), 1SHIZ(MMX),BNTR(MMX),CNTR(MMX),BRTR,QQOP(MMX) 2,XLOP(MMX),SHA(MMX),SHAZ(MMX),RR(MMX),SS(MMX),DAM(MMX), 3DAMZ(MMX),DEDY(MMX),DRZ(MMX),DSZ(MMX),ANTR(MMX) include "sodium.h" COMMON/WINDTRF/UN(KMX),VN(KMX),WN(KMX),UI(KMX),VI(KMX),WI(KMX) 1,WGAM(KMX),DWGAM(KMX),UIT,VUT C FLUXE=0. FLUXE=-5.E+6 M=KMX-MMX+1 BL(1)=XNASP(M) BR(1)=BRTR*AR(1)-FLUXE/DA(MMX)*SHA(MMX) RETURN END SUBROUTINE SOLB (NDIM, N, ML, MU, B, Y, IP) SAVE DIMENSION B(NDIM,1),Y(N),IP(N) C----------------------------------------------------------------------- C SOLUTION OF A*X = C GIVEN LU DECOMPOSITION OF A FROM DECB. C Y = RIGHT-HAND VECTOR C, OF LENGTH N, ON INPUT, C = SOLUTION VECTOR X ON OUTPUT. C ALL THE ARGUMENTS ARE INPUT ARGUMENTS. C THE OUTPUT ARGUMENT IS Y. C----------------------------------------------------------------------- IF (N .EQ. 1) GO TO 60 N1 = N - 1 LL = ML + MU + 1 IF (ML .EQ. 0) GO TO 32 DO 30 NR = 1,N1 IF (IP(NR) .EQ. NR) GO TO 10 J = IP(NR) XX = Y(NR) Y(NR) = Y(J) Y(J) = XX 10 KK = MIN0(N-NR,ML) DO 20 I = 1,KK 20 Y(NR+I) = Y(NR+I) + Y(NR)*B(NR,LL+I) 30 CONTINUE 32 LL = LL - 1 Y(N) = Y(N)*B(N,1) KK = 0 DO 50 NB = 1,N1 NR = N - NB IF (KK .NE. LL) KK = KK + 1 DP = 0. IF (LL .EQ. 0) GO TO 50 DO 40 I = 1,KK 40 DP = DP + B(NR,I+1)*Y(NR+I) 50 Y(NR) = (Y(NR) - DP)*B(NR,1) RETURN 60 Y(1) = Y(1)*B(1,1) RETURN C----------------------- END OF SUBROUTINE SOLB ------------------------ END SUBROUTINE DECB (NDIM, N, ML, MU, B, IP, IER) SAVE DIMENSION B(NDIM,1),IP(N) C----------------------------------------------------------------------- C LU DECOMPOSITION OF BAND MATRIX A.. L*U = P*A , WHERE P IS A C PERMUTATION MATRIX, L IS A UNIT LOWER TRIANGULAR MATRIX, C AND U IS AN UPPER TRIANGULAR MATRIX. C N = ORDER OF MATRIX. C B = N BY (2*ML+MU+1) ARRAY CONTAINING THE MATRIX A ON INPUT C AND ITS FACTORED FORM ON OUTPUT. C ON INPUT, B(I,K) (1.LE.I.LE.N) CONTAINS THE K-TH C DIAGONAL OF A, OR A(I,J) IS STORED IN B(I,J-I+ML+1). C ON OUTPUT, B CONTAINS THE L AND U FACTORS, WITH C U IN COLUMNS 1 TO ML+MU+1, AND L IN COLUMNS C ML+MU+2 TO 2*ML+MU+1. C ML,MU = WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, NOT C COUNTING THE MAIN DIAGONAL. TOTAL BANDWIDTH IS ML+MU+1. C NDIM = THE FIRST DIMENSION (COLUMN LENGTH) OF THE ARRAY B. C NDIM MUST BE .GE. N. C IP = ARRAY OF LENGTH N CONTAINING PIVOT INFORMATION. C IER = ERROR INDICATOR.. C = 0 IF NO ERROR, C = K IF THE K-TH PIVOT CHOSEN WAS ZERO (A IS SINGULAR). C THE INPUT ARGUMENTS ARE NDIM, N, ML, MU, B. C THE OUTPUT ARGUMENTS ARE B, IP, IER. C----------------------------------------------------------------------- IER = 0 IF (N .EQ. 1) GO TO 92 LL = ML + MU + 1 N1 = N - 1 IF (ML .EQ. 0) GO TO 32 DO 30 I = 1,ML II = MU + I K = ML + 1 - I DO 10 J = 1,II 10 B(I,J) = B(I,J+K) K = II + 1 DO 20 J = K,LL 20 B(I,J) = 0. 30 CONTINUE 32 LR = ML DO 90 NR = 1,N1 NP = NR + 1 IF (LR .NE. N) LR = LR + 1 MX = NR XM = ABS(B(NR,1)) IF (ML .EQ. 0) GO TO 42 DO 40 I = NP,LR IF (ABS(B(I,1)) .LE. XM) GO TO 40 MX = I XM = ABS(B(I,1)) 40 CONTINUE 42 IP(NR) = MX IF (MX .EQ. NR) GO TO 60 DO 50 I = 1,LL XX = B(NR,I) B(NR,I) = B(MX,I) 50 B(MX,I) = XX 60 XM = B(NR,1) IF (XM .EQ. 0.) GO TO 100 B(NR,1) = 1./XM IF (ML .EQ. 0) GO TO 90 XM = -B(NR,1) KK = MIN0(N-NR,LL-1) DO 80 I = NP,LR J = LL + I - NR XX = B(I,1)*XM B(NR,J) = XX DO 70 II = 1,KK 70 B(I,II) = B(I,II+1) + XX*B(NR,II+1) 80 B(I,LL) = 0. 90 CONTINUE 92 NR = N IF (B(N,1) .EQ. 0.) GO TO 100 B(N,1) = 1./B(N,1) RETURN 100 IER = NR RETURN C----------------------- END OF SUBROUTINE DECB ------------------------ END SUBROUTINE NAFUNC(NPDE,NZP2,ZN,TM,UM,UZN,AN,BN,CN) SAVE PARAMETER(KMX=96,MMX=96) include "blank.h" COMMON/TRAVLH/RK(50),ALP(20),RB(20),RKM(60) COMMON/COLFACT/COLFAC REAL ZN(NZP2),UM(NPDE,NZP2),AN(NPDE,NZP2),BN(NPDE,NZP2), 1CN(NPDE,NZP2),UZN(NPDE,NZP2) COMMON/SOLVT/TP(MMX),DA(MMX),SHI(MMX),DAZ(MMX), 1SHIZ(MMX),BNTR(MMX),CNTR(MMX),BRTR,QQOP(MMX) 2,XLOP(MMX),SHA(MMX),SHAZ(MMX),RR(MMX),SS(MMX),DAM(MMX), 3DAMZ(MMX),DEDY(MMX),DRZ(MMX),DSZ(MMX),ANTR(MMX) include "sodium.h" COMMON/WINDTRF/UN(KMX),VN(KMX),WN(KMX),UI(KMX),VI(KMX),WI(KMX) 1,WGAM(KMX),DWGAM(KMX),UIT,VUT MMX1=MMX+1 MMX2=MMX+2 M1=KMX-MMX DO 1 KK=2,MMX1 K=KK-1 M=M1+K AN(1,KK)=ANTR(K) BN(1,KK)=BNTR(K) CN(1,KK)=PNASP(K)+(CNTR(K)-XLNASP(K))*UM(1,KK) 1 CONTINUE AN(1,1)=AN(1,2) BN(1,1)=BN(1,2) CN(1,1)=CN(1,2) AN(1,MMX2)=AN(1,MMX1) BN(1,MMX2)=BN(1,MMX1) CN(1,MMX2)=CN(1,MMX1) RETURN END SUBROUTINE DIFFN(NPDE,NZP2,NEQ,T,ZN,U,UT,UT1,DM,UN,AN,BN,CN,AL, 1 AR,BL,BR,MBLFT,MBRHT,UZN,FUNC,BNBDY) SAVE C . C THIS ROUTINE COMPUTES THE TIME DERIVATIVE GIVEN THE SPATIAL C DERIVATIVES. C . REAL ZN(NZP2),UN(NPDE,NZP2),U(NEQ),UT(NEQ),UT1(NZP2), 1 BN(NPDE,NZP2),CN(NPDE,NZP2),AL(NPDE),AR(NPDE),DM(NZP2,3), 2 AN(NPDE,NZP2),BL(NPDE),BR(NPDE),UZN(NPDE,NZP2) INTEGER MBLFT(NPDE),MBRHT(NPDE) COMMON /PDVAR/ NPDEC,NZC,NBUG,ITYPE COMMON /SYSPAR/ UROUND,MIN,MOUT EXTERNAL FUNC,BNBDY C . NZ=NZP2-2 NZP1=NZP2-1 C . C MOVE THE SOLUTION IN THE U ARRAY TO THE ARRAY UN AND SET THE C BOUNDARY VALUES FROM THE BOUNDARY CONDITIONS C . CALL MVTOUN(T,NPDE,NZP2,NEQ,U,MBLFT,MBRHT,BL,BR,UN,ZN,BNBDY) C . C COMPUTE THE TIME DERIVATIVE AND STORE TEMPORARILY IN THE UT1 ARRAY C . C . C DETERMINE THE NUMBER OF MIXED BOUNDARY VARIABLES AT THE LEFT C (NQLFT) AND AT THE RIGHT (NQRHT) C . NQLFT=0 NQRHT=0 DO 50 M=1,NPDE IF(MBLFT(M) .NE. 1)NQLFT=NQLFT+1 IF(MBRHT(M) .NE. 1)NQRHT=NQRHT+1 50 CONTINUE C . C COMPUTE THE DERIVATIVE UZN BY 2-ND ORDER DIFFERENCES C . DO 110 I=2,NZP1 DM(I,1)=1./((ZN(I)-ZN(I-1))*(ZN(I+1)-ZN(I-1))) DM(I,2)=1./((ZN(I)-ZN(I-1))*(ZN(I+1)-ZN(I))) DM(I,3)=1./((ZN(I+1)-ZN(I))*(ZN(I+1)-ZN(I-1))) 110 CONTINUE DO 120 M=1,NPDE DO 120 I=2,NZP1 UZN(M,I)= -(ZN(I+1)-ZN(I))*UN(M,I-1)*DM(I,1) 1 +(ZN(I+1)-2.*ZN(I)+ZN(I-1))*UN(M,I)*DM(I,2) 2 +(ZN(I)-ZN(I-1))*UN(M,I+1)*DM(I,3) 120 CONTINUE C C COMPUTE THE TIME DERIVATVE AT THE BOUNDARY IF THERE IS A TYPE C THREE BOUNDARY CONDITION (MBLFT=3 OR MBRHT=3) C DO 130 M=1,NPDE IF(MBLFT(M) .EQ. 3)GOTO 140 IF(MBRHT(M) .EQ. 3)GOTO 140 130 CONTINUE GOTO 150 140 CALL BNBDY(NPDE,T,ZN(2),ZN(NZP1),UN(1,2),UN(1,NZP1), 1 BL,BR) 150 CONTINUE C . C EVALUATE COEFFICIENT FUNCTIONS C . CALL FUNC(NPDE,NZP2,ZN,T,UN,UZN,AN,BN,CN) C . C DEBUG PRINTOUT NBUG=8 C . IF(NBUG/8-2*(NBUG/16) .NE. 1)GOTO 190 WRITE (MOUT,180) T 180 FORMAT(//3X,24HNBUG=2 PRINTING IN DIFFN,6H AT T=, 1 E12.5) DO 185 M=1,NPDE WRITE (MOUT,181) M 181 FORMAT(/3X,2HM=,I3/5X,1HI,9X,2HZN,9X,2HUN,8X,3HUZN, 1 9X,2HAN,9X,2HBN,9X,2HCN) DO 185 I=1,NZP2 WRITE (MOUT,182) I,ZN(I),UN(M,I),UZN(M,I),AN(M,I),BN(M,I), 1 CN(M,I) 182 FORMAT(1X,I3,2X,6E12.4) 185 CONTINUE 190 CONTINUE C . KTL=0 KTR=NPDE*(NZ-2) + NQLFT DO 280 M=1,NPDE IF(ITYPE .EQ. 2)GOTO 222 C . C CONSERVATION FORM OF THE SYSTEM C . DO 220 I=2,NZP1 T1=.5*(AN(M,I+1)+AN(M,I))*(UN(M,I+1)-UN(M,I)) T2=T1/(ZN(I+1)-ZN(I)) T3=.5*(AN(M,I)+AN(M,I-1))*(UN(M,I)-UN(M,I-1)) T4=T3/(ZN(I)-ZN(I-1)) T5=.25*(BN(M,I+1)+BN(M,I))*(UN(M,I+1)+UN(M,I)) T6=.25*(BN(M,I)+BN(M,I-1))*(UN(M,I)+UN(M,I-1)) T7=(T2-T4+T5-T6)/(.5*(ZN(I+1)-ZN(I-1))) T8=T7+CN(M,I) UT1(I)=T8 220 CONTINUE GOTO 228 C . C NON CONSERVATION FORM OF THE EQUATIONS C . 222 DO 224 I=2,NZP1 UT1(I)=AN(M,I)*2.*(UN(M,I-1)*DM(I,1) - UN(M,I)*DM(I,2) 1 +UN(M,I+1)*DM(I,3)) 2 +BN(M,I)*(-(ZN(I+1)-ZN(I))*UN(M,I-1)*DM(I,1) 3 +(ZN(I+1)-2.*ZN(I)+ZN(I-1))*UN(M,I)*DM(I,2) 4 +(ZN(I)-ZN(I-1))*UN(M,I+1)*DM(I,3) ) 5 + CN(M,I) 224 CONTINUE 228 CONTINUE C C CORRECT TIME DERIVATIVE AT THE BOUNDARY C IF(MBLFT(M) .EQ. 3)UT1(2)=BL(M) IF(MBRHT(M) .EQ. 3)UT1(NZP1)=BR(M) C . C STORE THE TIME DERIVATIVE IN THE ARRAY UT FOR USE BY GEARB C . IF(MBLFT(M) .EQ. 1)GOTO 230 KTL=KTL+1 UT(KTL)=UT1(2) 230 IF(MBRHT(M) .EQ. 1)GOTO 240 KTR=KTR+1 UT(KTR)=UT1(NZP1) 240 KT=NQLFT+M DO 245 I=3,NZ UT(KT)=UT1(I) 245 KT=KT+NPDE C . C DEBUG PRINTOUT NBUG=8 C . IF(NBUG/8-2*(NBUG/16) .NE. 1)GOTO 280 WRITE (MOUT,270) T,M 270 FORMAT(/3X,24HNBUG=4 PRINTING IN DIFFN,6H AT T=, 1 E12.5,3X,2HM=,I3/10X,12HUT1(2..NZP1)) WRITE (MOUT,272) (UT1(I),I=2,NZP1) 272 FORMAT(1X,6E12.4) 280 CONTINUE C . RETURN END SUBROUTINE DIFFUN(NEQ,T,U,UT,WKCM,IWCM,NWKCM,NIWCM, 1 FUNC,BNBDY) SAVE C . C THE ROUTINE REQUIRED BY GEARB TO EVALUTE THE TIME DERIVATIVE C . REAL U(NEQ),UT(NEQ),WKCM(NWKCM) INTEGER IWCM(NIWCM) COMMON /PDVAR/ NPDE,NZ,NBUG,ITYPE COMMON /SYSPAR/ UROUND,MIN,MOUT EXTERNAL FUNC,BNBDY C . NZP2=NZ+2 LML=1 LMR=LML+NPDE LZN=1 LUN=LZN+NZP2 LU=LUN+NPDE*NZP2 LUT=LU+NPDE*NZP2 LUZN=LUT+NPDE*NZP2 LUT1=LUZN+NPDE*NZP2 LDM=LUT1+NZP2 LAN=LDM+3*NZP2 LBN=LAN+NPDE*NZP2 LCN=LBN+NPDE*NZP2 LAL=LCN+NPDE*NZP2 LAR=LAL+NPDE LBL=LAR+NPDE LBR=LBL+NPDE CALL DIFFN(NPDE,NZP2,NEQ,T,WKCM(LZN),U,UT,WKCM(LUT1),WKCM(LDM), 1 WKCM(LUN),WKCM(LAN),WKCM(LBN),WKCM(LCN),WKCM(LAL),WKCM(LAR), 2 WKCM(LBL),WKCM(LBR),IWCM(LML),IWCM(LMR),WKCM(LUZN),FUNC,BNBDY) C . C DEBUG PRINTOUT WHEN NBUG=4 C . IF(NBUG/4-2*(NBUG/8) .NE. 1)GOTO 75 WRITE (MOUT,65) T 65 FORMAT(//3X,25HNBUG=1 PRINTING IN DIFFUN,6H AT T=, 1 E12.5/5X,1HK,8X,4HU(I),7X,5HUT(I)) DO 70 K=1,NEQ WRITE (MOUT,67) K,U(K),UT(K) 67 FORMAT(3X,I3,3X,2E13.5) 70 CONTINUE 75 CONTINUE RETURN END SUBROUTINE COSET (METH, NQ, EL, TQ, MAXDER) SAVE C----------------------------------------------------------------------- C COSET IS CALLED BY THE INTEGRATOR AND SETS COEFFICIENTS USED THERE. C THE VECTOR EL, OF LENGTH NQ + 1, DETERMINES THE BASIC METHOD. C THE VECTOR TQ, OF LENGTH 4, IS INVOLVED IN ADJUSTING THE STEP SIZE C IN RELATION TO TRUNCATION ERROR. ITS VALUES ARE GIVEN BY THE C PERTST ARRAY. C THE VECTORS EL AND TQ DEPEND ON METH AND NQ. C COSET ALSO SETS MAXDER, THE MAXIMUM ORDER OF THE METHOD AVAILABLE. C CURRENTLY IT IS 12 FOR THE ADAMS METHODS AND 5 FOR THE BDF METHODS. C LMAX = MAXDER + 1 IS THE NUMBER OF COLUMNS IN THE Y ARRAY. C THE MAXIMUM ORDER USED MAY BE REDUCED SIMPLY BY DECREASING THE C NUMBERS IN STATEMENTS 1 AND/OR 2 BELOW. C . C THE COEFFICIENTS IN PERTST NEED BE GIVEN TO ONLY ABOUT C ONE PERCENT ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW C IS.. COEFFICIENTS FOR ORDER NQ - 1, COEFFICIENTS FOR ORDER NQ, C COEFFICIENTS FOR ORDER NQ + 1. WITHIN EACH GROUP ARE THE C COEFFICIENTS FOR THE ADAMS METHODS, FOLLOWED BY THOSE FOR THE C BDF METHODS. C----------------------------------------------------------------------- DIMENSION PERTST(12,2,3),EL(13),TQ(4) DATA PERTST / 1.,1.,2.,1.,.3158,.07407,.01391,.002182, 1 .0002945,.00003492,.000003692,.0000003524, 2 1.,1.,.5,.1667,.04167,1.,1.,1.,1.,1.,1.,1., 3 2.,12.,24.,37.89,53.33,70.08,87.97,106.9, 4 126.7,147.4,168.8,191.0, 5 2.0,4.5,7.333,10.42,13.7,1.,1.,1.,1.,1.,1.,1., 6 12.0,24.0,37.89,53.33,70.08,87.97,106.9, 7 126.7,147.4,168.8,191.0,1., 8 3.0,6.0,9.167,12.5,1.,1.,1.,1.,1.,1.,1.,1. / C . GO TO (1,2),METH 1 MAXDER = 12 GO TO (101,102,103,104,105,106,107,108,109,110,111,112),NQ 2 MAXDER = 5 GO TO (201,202,203,204,205),NQ C----------------------------------------------------------------------- C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY. C FOR A GIVEN ORDER NQ, THEY CAN BE CALCULATED BY USE OF THE C GENERATING POLYNOMIAL L(T), WHOSE COEFFICIENTS ARE EL(I).. C L(T) = EL(1) + EL(2)*T + ... + EL(NQ+1)*T**NQ. C FOR THE IMPLICIT ADAMS METHODS, L(T) IS GIVEN BY C DL/DT = (T+1)*(T+2)* ... *(T+NQ-1)/K, L(-1) = 0, C WHERE K = FACTORIAL(NQ-1). C FOR THE BDF METHODS, C L(T) = (T+1)*(T+2)* ... *(T+NQ)/K, C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ). C . C THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS.. C IMPLICIT ADAMS METHODS OF ORDERS 1 TO 12, C BDF METHODS OF ORDERS 1 TO 5. C----------------------------------------------------------------------- 101 EL(1) = 1.0 GO TO 900 102 EL(1) = 0.5 EL(3) = 0.5 GO TO 900 103 EL(1) = 4.1666666666667E-01 EL(3) = 0.75 EL(4) = 1.6666666666667E-01 GO TO 900 104 EL(1) = 0.375 EL(3) = 9.1666666666667E-01 EL(4) = 3.3333333333333E-01 EL(5) = 4.1666666666667E-02 GO TO 900 105 EL(1) = 3.4861111111111E-01 EL(3) = 1.0416666666667 EL(4) = 4.8611111111111E-01 EL(5) = 1.0416666666667E-01 EL(6) = 8.3333333333333E-03 GO TO 900 106 EL(1) = 3.2986111111111E-01 EL(3) = 1.1416666666667 EL(4) = 0.625 EL(5) = 1.7708333333333E-01 EL(6) = 0.025 EL(7) = 1.3888888888889E-03 GO TO 900 107 EL(1) = 3.1559193121693E-01 EL(3) = 1.225 EL(4) = 7.5185185185185E-01 EL(5) = 2.5520833333333E-01 EL(6) = 4.8611111111111E-02 EL(7) = 4.8611111111111E-03 EL(8) = 1.9841269841270E-04 GO TO 900 108 EL(1) = 3.0422453703704E-01 EL(3) = 1.2964285714286 EL(4) = 8.6851851851852E-01 EL(5) = 3.3576388888889E-01 EL(6) = 7.7777777777778E-02 EL(7) = 1.0648148148148E-02 EL(8) = 7.9365079365079E-04 EL(9) = 2.4801587301587E-05 GO TO 900 109 EL(1) = 2.9486800044092E-01 EL(3) = 1.3589285714286 EL(4) = 9.7655423280423E-01 EL(5) = 0.4171875 EL(6) = 1.1135416666667E-01 EL(7) = 0.01875 EL(8) = 1.9345238095238E-03 EL(9) = 1.1160714285714E-04 EL(10)= 2.7557319223986E-06 GO TO 900 110 EL(1) = 2.8697544642857E-01 EL(3) = 1.4144841269841 EL(4) = 1.0772156084656 EL(5) = 4.9856701940035E-01 EL(6) = 0.1484375 EL(7) = 2.9060570987654E-02 EL(8) = 3.7202380952381E-03 EL(9) = 2.9968584656085E-04 EL(10)= 1.3778659611993E-05 EL(11)= 2.7557319223986E-07 GO TO 900 111 EL(1) = 2.8018959644394E-01 EL(3) = 1.4644841269841 EL(4) = 1.1715145502646 EL(5) = 5.7935819003527E-01 EL(6) = 1.8832286155203E-01 EL(7) = 4.1430362654321E-02 EL(8) = 6.2111441798942E-03 EL(9) = 6.2520667989418E-04 EL(10)= 4.0417401528513E-05 EL(11)= 1.5156525573192E-06 EL(12)= 2.5052108385442E-08 GO TO 900 112 EL(1) = 2.7426554003160E-01 EL(3) = 1.5099386724387 EL(4) = 1.2602711640212 EL(5) = 6.5923418209877E-01 EL(6) = 2.3045800264550E-01 EL(7) = 5.5697246105232E-02 EL(8) = 9.4394841269841E-03 EL(9) = 1.1192749669312E-03 EL(10)= 9.0939153439153E-05 EL(11)= 4.8225308641975E-06 EL(12)= 1.5031265031265E-07 EL(13)= 2.0876756987868E-09 GO TO 900 201 EL(1) = 1.0 GO TO 900 202 EL(1) = 6.6666666666667E-01 EL(3) = 3.3333333333333E-01 GO TO 900 203 EL(1) = 5.4545454545455E-01 EL(3) = EL(1) EL(4) = 9.0909090909091E-02 GO TO 900 204 EL(1) = 0.48 EL(3) = 0.7 EL(4) = 0.2 EL(5) = 0.02 GO TO 900 205 EL(1) = 4.3795620437956E-01 EL(3) = 8.2116788321168E-01 EL(4) = 3.1021897810219E-01 EL(5) = 5.4744525547445E-02 EL(6) = 3.6496350364964E-03 C . 900 DO 910 K = 1,3 910 TQ(K) = PERTST(NQ,METH,K) TQ(4) = .5*TQ(2)/FLOAT(NQ+2) RETURN C----------------------- END OF SUBROUTINE COSET ----------------------- END