SUBROUTINE INITAL C... C... H5 - DYNAMIC GRAETZ PROBLEMS WITH CONSTANT WALL TEMPERATURE C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... SELECT APPROXIMATION OF SPATIAL DERIVATIVES (SEE THE COMMENTS C... IN SUBROUTINE DERV FOR AN EXPLANATION OF THE EIGHT CASES) NCASE=8 C... C... PROBLEM PARAMETERS PE=60.0D0 R0= 1.0D0 IF(NORUN.EQ.1)ZL=30.0D0 IF(NORUN.EQ.2)ZL=90.0D0 C... C... RADIAL GRID, VELOCITY PROFILE DR=R0/DFLOAT(NR-1) DO 1 IR=1,NR R(IR)=DFLOAT(IR-1)*DR V(IR)=2.0D0*(1.0D0-R(IR)**2) 1 CONTINUE DRS=DR**2 C... C... AXIAL GRID DZ=ZL/DFLOAT(NZ-1) DO 2 JZ=1,NZ Z(JZ)=DFLOAT(JZ-1)*DZ 2 CONTINUE DZS=DZ**2 C... C... INITIAL CONDITION DO 3 IR=1,NR DO 3 JZ=1,NZ U(IR,JZ)=0.0D0 3 CONTINUE C... C... GIVE THE TIME DERIVATIVES A VALUE WHICH CAN THEN BE USED TO C... TEST FOR UNDEFINED DERIVATIVES BY SETTING THEM TO PI, THEN C... CALLING DERV, AND CHECKING IN SUBROUTINE PRINT IF THEY HAVE C... ALL BEEN SET TO VALUES OTHER THAN PI. THIS IS A SIMPLE PRO- C... CEDURE TO ENSURE THAT ALL NR X NZ DERIVATIVES HAVE BEEN SET C... IN DERV PI=3.1415927D0 DO 4 IR=1,NR DO 4 JZ=1,NZ UT(IR,JZ)=PI 4 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/I/ IP, NCASE C... C... APPROXIMATIONS OF SPATIAL DERIVATIVES C... C... NCASE = 1 - EXPLICIT FINITE DIFFERENCES (FDM) WITH UZZ, C... UZ(R,ZL,T) = 0 IF(NCASE.EQ.1)THEN CALL DERV1 C... C... NCASE = 2 - EXPLICIT FINITE DIFFERENCES WITH UZZ, C... UT(R,ZL,T) = -PE*V(R)*UZ(R,ZL,T) ELSE + IF(NCASE.EQ.2)THEN CALL DERV2 C... C... NCASE = 3 - EXPLICIT FINITE DIFFERENCES WITH UZZ, C... UT(R,ZL,T) = -PE*V(R)*UZ(R,ZL,T) + C... URR(R,ZL,T) + (1/R)*UR(R,ZL,T) ELSE + IF(NCASE.EQ.3)THEN CALL DERV3 C... C... NCASE = 4 - EXPLICIT FINITE DIFFERENCES WITHOUT UZZ ELSE + IF(NCASE.EQ.4)THEN CALL DERV4 C... C... NCASE = 5 - FINITE DIFFERENCES FROM DSS034 WITH UZZ, C... UZ(R,ZL,T) = 0 ELSE + IF(NCASE.EQ.5)THEN CALL DERV5 C... C... NCASE = 6 - FINITE DIFFERENCES FROM DSS034 WITH UZZ, C... UT(R,ZL,T) = -PE*V(R)*UZ(R,ZL,T) ELSE + IF(NCASE.EQ.6)THEN CALL DERV6 C... C... NCASE = 7 - FINITE DIFFERENCES FROM DSS034 WITH UZZ, C... UT(R,ZL,T) = -PE*V(R)*UZ(R,ZL,T) + C... URR(R,ZL,T) + (1/R)*UR(R,ZL,T) ELSE + IF(NCASE.EQ.7)THEN CALL DERV7 C... C... NCASE = 8 - FINITE DIFFERENCES FROM DSS034 WITHOUT UZZ ELSE + IF(NCASE.EQ.8)THEN CALL DERV8 END IF RETURN END SUBROUTINE DERV1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=0.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ U(NR,JZ)=1.0D0 UT(NR,JZ)=0.0D0 2 CONTINUE C... C... R = 0, Z NE 0, ZL DO 3 JZ=2,NZ-1 UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ + +4.0D0*(U(2,JZ)-U(1,JZ))/DRS + +(U(1,JZ+1)-2.0D0*U(1,JZ)+U(1,JZ-1))/DZS 3 CONTINUE C... C... R = 0, Z = ZL UT(1,NZ)=-PE*V(1)*0.0D0 + +4.0D0*(U(2,NZ) -U(1,NZ))/DRS + +2.0D0*(U(1,NZ-1)-U(1,NZ))/DZS C... C... Z = ZL, R NE 0, R0 DO 4 IR=2,NR-1 UT(IR,NZ)=-PE*V(IR)*0.0D0 + +(U(IR+1,NZ)-2.0D0*U(IR,NZ)+U(IR-1,NZ))/DRS + +(1.0D0/R(IR))*(U(IR+1,NZ)-U(IR-1,NZ))/(2.0D0*DR) + +2.0D0*(U(IR,NZ-1)-U(IR,NZ))/DZS 4 CONTINUE C... C... R NE 0, R0, Z NE 0, ZL DO 5 JZ=2,NZ-1 DO 5 IR=2,NR-1 UT(IR,JZ)=-PE*V(IR)*(U(IR,JZ)-U(IR,JZ-1))/DZ + +(U(IR+1,JZ)-2.0D0*U(IR,JZ)+U(IR-1,JZ))/DRS + +(1.0D0/R(IR))*(U(IR+1,JZ)-U(IR-1,JZ))/(2.0D0*DR) + +(U(IR,JZ+1)-2.0D0*U(IR,JZ)+U(IR,JZ-1))/DZS 5 CONTINUE RETURN END SUBROUTINE DERV2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=0.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ U(NR,JZ)=1.0D0 UT(NR,JZ)=0.0D0 2 CONTINUE C... C... R = 0, Z NE 0, ZL DO 3 JZ=2,NZ-1 UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ + +4.0D0*(U(2,JZ)-U(1,JZ))/DRS + +(U(1,JZ+1)-2.0D0*U(1,JZ)+U(1,JZ-1))/DZS 3 CONTINUE C... C... Z = ZL, R NE R0 DO 4 IR=1,NR-1 UT(IR,NZ)=-PE*V(IR)*(U(IR,NZ)-U(IR,NZ-1))/DZ 4 CONTINUE C... C... R NE 0, R0, Z NE 0, ZL DO 5 JZ=2,NZ-1 DO 5 IR=2,NR-1 UT(IR,JZ)=-PE*V(IR)*(U(IR,JZ)-U(IR,JZ-1))/DZ + +(U(IR+1,JZ)-2.0D0*U(IR,JZ)+U(IR-1,JZ))/DRS + +(1.0D0/R(IR))*(U(IR+1,JZ)-U(IR-1,JZ))/(2.0D0*DR) + +(U(IR,JZ+1)-2.0D0*U(IR,JZ)+U(IR,JZ-1))/DZS 5 CONTINUE RETURN END SUBROUTINE DERV3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=0.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ U(NR,JZ)=1.0D0 UT(NR,JZ)=0.0D0 2 CONTINUE C... C... R = 0, Z NE 0, ZL DO 3 JZ=2,NZ-1 UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ + +4.0D0*(U(2,JZ)-U(1,JZ))/DRS + +(U(1,JZ+1)-2.0D0*U(1,JZ)+U(1,JZ-1))/DZS 3 CONTINUE C... C... R = 0, Z = ZL UT(1,NZ)=-PE*V(1)*(U(1,NZ)-U(1,NZ-1))/DZ + +4.0D0*(U(2,NZ) -U(1,NZ))/DRS C... C... Z = ZL, R NE 0, R0 DO 4 IR=2,NR-1 UT(IR,NZ)=-PE*V(IR)*(U(IR,NZ)-U(IR,NZ-1))/DZ + +(U(IR+1,NZ)-2.0D0*U(IR,NZ)+U(IR-1,NZ))/DRS + +(1.0D0/R(IR))*(U(IR+1,NZ)-U(IR-1,NZ))/(2.0D0*DR) 4 CONTINUE C... C... R NE 0, R0, Z NE 0, ZL DO 5 JZ=2,NZ-1 DO 5 IR=2,NR-1 UT(IR,JZ)=-PE*V(IR)*(U(IR,JZ)-U(IR,JZ-1))/DZ + +(U(IR+1,JZ)-2.0D0*U(IR,JZ)+U(IR-1,JZ))/DRS + +(1.0D0/R(IR))*(U(IR+1,JZ)-U(IR-1,JZ))/(2.0D0*DR) + +(U(IR,JZ+1)-2.0D0*U(IR,JZ)+U(IR,JZ-1))/DZS 5 CONTINUE RETURN END SUBROUTINE DERV4 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=0.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ U(NR,JZ)=1.0D0 UT(NR,JZ)=0.0D0 2 CONTINUE C... C... R = 0, Z NE 0 DO 3 JZ=2,NZ UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ + +4.0D0*(U(2,JZ)-U(1,JZ))/DRS 3 CONTINUE C... C... R NE 0, R0, Z NE 0 DO 5 JZ=2,NZ DO 5 IR=2,NR-1 UT(IR,JZ)=-PE*V(IR)*(U(IR,JZ)-U(IR,JZ-1))/DZ + +(U(IR+1,JZ)-2.0D0*U(IR,JZ)+U(IR-1,JZ))/DRS + +(1.0D0/R(IR))*(U(IR+1,JZ)-U(IR-1,JZ))/(2.0D0*DR) 5 CONTINUE RETURN END SUBROUTINE DERV5 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... ARRAYS FOR UR, URR, UZ, UZZ DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ), UZZ(NR,NZ) C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=0.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ U(NR,JZ)=1.0D0 UT(NR,JZ)=0.0D0 2 CONTINUE C... C... UR CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0) C... C... R = 0 DO 3 JZ=1,NZ UR(1,JZ)=0.0D0 3 CONTINUE C... C... URR CALL DSS034(0.0D0,R0,NR,NZ,1,UR,URR,0.0D0) C... C... UZ (FOR CONDUCTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,0.0D0) C... C... Z = ZL DO 4 IR=1,NR UZ(IR,NZ)=0.0D0 4 CONTINUE C... C... UZZ (FOR CONDUCTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,UZ,UZZ,0.0D0) C... C... UZ (FOR CONVECTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0) C... C... Z = ZL DO 5 IR=1,NR UZ(IR,NZ)=0.0D0 5 CONTINUE C... C... R NE R0, Z NE 0 DO 6 JZ=2,NZ DO 6 IR=1,NR-1 C... C... R = 0 IF(IR.EQ.1)THEN UT(1,JZ)=-PE*V(1)*UZ(1,JZ) + +2.0D0*URR(1,JZ) + +UZZ(1,JZ) C... C... R NE 0 ELSE + IF(IR.NE.1)THEN UT(IR,JZ)=-PE*V(IR)*UZ(IR,JZ) + +URR(IR,JZ)+(1.0D0/R(IR))*UR(IR,JZ) + +UZZ(IR,JZ) END IF 6 CONTINUE RETURN END SUBROUTINE DERV6 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... ARRAYS FOR UR, URR, UZ, UZZ DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ), UZZ(NR,NZ) C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=0.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ U(NR,JZ)=1.0D0 UT(NR,JZ)=0.0D0 2 CONTINUE C... C... UR CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0) C... C... R = 0 DO 3 JZ=1,NZ UR(1,JZ)=0.0D0 3 CONTINUE C... C... URR CALL DSS034(0.0D0,R0,NR,NZ,1,UR,URR,0.0D0) C... C... UZ (FOR CONDUCTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,0.0D0) C... C... UZZ (FOR CONDUCTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,UZ,UZZ,0.0D0) C... C... UZ (FOR CONVECTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0) C... C... R NE R0, Z NE 0, ZL DO 4 JZ=2,NZ-1 DO 4 IR=1,NR-1 C... C... R = 0 IF(IR.EQ.1)THEN UT(1,JZ)=-PE*V(1)*UZ(1,JZ) + +2.0D0*URR(1,JZ) + +UZZ(1,JZ) C... C... R NE 0 ELSE + IF(IR.NE.1)THEN UT(IR,JZ)=-PE*V(IR)*UZ(IR,JZ) + +URR(IR,JZ)+(1.0D0/R(IR))*UR(IR,JZ) + +UZZ(IR,JZ) END IF 4 CONTINUE C... C... R NE R0, Z = ZL DO 5 IR=1,NR-1 UT(IR,NZ)=-PE*V(IR)*UZ(IR,NZ) 5 CONTINUE RETURN END SUBROUTINE DERV7 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... ARRAYS FOR UR, URR, UZ, UZZ DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ), UZZ(NR,NZ) C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=0.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ U(NR,JZ)=1.0D0 UT(NR,JZ)=0.0D0 2 CONTINUE C... C... UR CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0) C... C... R = 0 DO 3 JZ=1,NZ UR(1,JZ)=0.0D0 3 CONTINUE C... C... URR CALL DSS034(0.0D0,R0,NR,NZ,1,UR,URR,0.0D0) C... C... UZ (FOR CONDUCTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,0.0D0) C... C... UZZ (FOR CONDUCTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,UZ,UZZ,0.0D0) C... C... UZ (FOR CONVECTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0) C... C... R NE R0, Z NE 0, ZL DO 4 JZ=2,NZ-1 DO 4 IR=1,NR-1 C... C... R = 0 IF(IR.EQ.1)THEN UT(1,JZ)=-PE*V(1)*UZ(1,JZ) + +2.0D0*URR(1,JZ) + +UZZ(1,JZ) C... C... R NE 0 ELSE + IF(IR.NE.1)THEN UT(IR,JZ)=-PE*V(IR)*UZ(IR,JZ) + +URR(IR,JZ)+(1.0D0/R(IR))*UR(IR,JZ) + +UZZ(IR,JZ) END IF 4 CONTINUE C... C... R NE R0, Z = ZL DO 5 IR=1,NR-1 C... C... R = 0 IF(IR.EQ.1)THEN UT(1,NZ)=-PE*V(1)*UZ(1,NZ) + +2.0D0*URR(1,NZ) C... C... R NE 0 ELSE + IF(IR.NE.1)THEN UT(IR,NZ)=-PE*V(IR)*UZ(IR,NZ) + +URR(IR,NZ)+(1.0D0/R(IR))*UR(IR,NZ) END IF 5 CONTINUE RETURN END SUBROUTINE DERV8 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... ARRAYS FOR UR, URR, UZ, UZZ DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ), UZZ(NR,NZ) C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=0.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ U(NR,JZ)=1.0D0 UT(NR,JZ)=0.0D0 2 CONTINUE C... C... UR CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0) C... C... R = 0 DO 3 JZ=1,NZ UR(1,JZ)=0.0D0 3 CONTINUE C... C... URR CALL DSS034(0.0D0,R0,NR,NZ,1,UR,URR,0.0D0) C... C... UZ (FOR CONVECTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0) C... C... R NE R0 DO 4 JZ=2,NZ DO 4 IR=1,NR-1 C... C... R = 0 IF(IR.EQ.1)THEN UT(1,JZ)=-PE*V(1)*UZ(1,JZ) + +2.0D0*URR(1,JZ) C... C... R NE 0 ELSE + IF(IR.NE.1)THEN UT(IR,JZ)=-PE*V(IR)*UZ(IR,JZ) + +URR(IR,JZ)+(1.0D0/R(IR))*UR(IR,JZ) END IF 4 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... MONITOR SOLUTION WRITE(*,*)NCASE,NORUN,T C... C... PRINT THE INITIAL CONDITION AND INITIAL TEMPORAL DERIVATIVE IF(IP.EQ.0)THEN WRITE(NO,10)T,((U(IR,JZ),IR=1,NR),JZ=1,NZ) 10 FORMAT(' t = ',F6.2,//,' u(r,z,t)',/,(11F6.3)) WRITE(NO,11)(( UT(IR,JZ),IR=1,NR),JZ=1,NZ) 11 FORMAT(//,' ut(r,z,t)',/,(11F6.1)) C... C... CHECK FOR UNDEFINED DERIVATIVES DO 12 IR=1,NR DO 12 JZ=1,NZ IF((UT(IR,JZ).GT.3.14D0).AND.(UT(IR,JZ).LT.3.15D0))THEN WRITE (NO,13)IR,JZ,U(IR,JZ),UT(IR,JZ) 13 FORMAT(' IR = ',I3,' JZ = ',I3, + ' U(IR,JZ) = ',F7.3,' UT(IR,JZ) = ',F7.1) END IF 12 CONTINUE END IF C... C... PRINT THE NUMERICAL SOLUTION EVERY FIFTH CALL TO PRINT IP=IP+1 IF(((IP-1)/5*5).EQ.(IP-1))THEN WRITE(NO,1)T 1 FORMAT(//,' t = ',F6.2,//,5X,'u(r,z,t)',//, + 4X,'z',4X,'r = 0.0',3X,'r = 0.2',3X,'r = 0.4', + 3X,'r = 0.6',3X,'r = 0.8',3X,'r = 1.0') DO 3 JZ=1,NZ NRI=2 WRITE(NO,2)Z(JZ),(U(IR,JZ),IR=1,NR,NRI) 2 FORMAT(F5.1,6F10.4) 3 CONTINUE END IF C... C... OPEN FILE FOR MATLAB PLOTTING IF((NORUN.EQ.1).AND.(T.LT.0.001D0))THEN OPEN(1,FILE='h5a.out') OPEN(2,FILE='h5b.out') END IF C... C... WRITE NUMERICAL SOLUTION AT R = 0 FOR MATLAB PLOTTING DO 4 JZ=1,NZ IF(NORUN.EQ.1)WRITE(1,5)Z(JZ),U(1,JZ) IF(NORUN.EQ.2)WRITE(2,5)Z(JZ),U(1,JZ) 5 FORMAT(2F10.4) 4 CONTINUE RETURN END