C ************************** C * PDE2D 9.2 MAIN PROGRAM * C ************************** C *** 3D PROBLEM SOLVED *** C############################################################################## C Is double precision mode to be used? Double precision is recommended # C on 32-bit computers. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If double precision mode is used, variables and functions assigned +# C + names beginning with a letter in the range A-H or O-Z will be DOUBLE +# C + PRECISION, and you should use double precision constants and FORTRAN +# C + expressions throughout; otherwise such variables and functions will +# C + be of type REAL. In either case, variables and functions assigned +# C + names beginning with I,J,K,L,M or N will be of INTEGER type. +# C + +# C + It is possible to convert a single precision PDE2D program to double +# C + precision after it has been created, using an editor. Just change +# C + all occurrences of "real" to "double precision" +# C + " tdp" to "dtdp" (note leading blank) +# C + Any user-written code or routines must be converted "by hand", of +# C + course. To convert from double to single, reverse the changes. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## implicit double precision (a-h,o-z) parameter (neqnmx= 99) C############################################################################## C NP1GRID = number of P1-grid lines # C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- PARAMETER (NP1GRID = 11) C############################################################################## C NP2GRID = number of P2-grid lines # C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- PARAMETER (NP2GRID = 11) C############################################################################## C NP3GRID = number of P3-grid lines # C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- PARAMETER (NP3GRID = 7) C############################################################################## C How many differential equations (NEQN) are there in your problem? # C############################################################################## PARAMETER (NEQN = 1) C DIMENSIONS OF WORK ARRAYS C SET TO 1 FOR AUTOMATIC ALLOCATION PARAMETER (IRWK8Z= 1) PARAMETER (IIWK8Z= 1) PARAMETER (NXP8Z=41,NYP8Z=41,NZP8Z=41,KDEG8Z=1) C############################################################################## C The solution is normally saved on an NP1+1 by NP2+1 by NP3+1 # C rectangular grid of points, # C P1 = P1A + I*(P1B-P1A)/NP1, I = 0,...,NP1 # C P2 = P2A + J*(P2B-P2A)/NP2, J = 0,...,NP2 # C P3 = P3A + K*(P3B-P3A)/NP3, K = 0,...,NP3 # C Enter values for NP1, NP2 and NP3. Suggested values: NP1=NP2=NP3=16. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you want to save the solution at an arbitrary user-specified set +# C + of points, set NP2=NP3=0 and NP1+1=number of points. In this case +# C + you can request tabular output, but no plots can be made. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## PARAMETER (NP1 = 40) PARAMETER (NP2 = 40) PARAMETER (NP3 = 40) PARAMETER (NSAVE = 1) common/parm8z/ pi,D dimension p1grid(np1grid),p2grid(np2grid),p3grid(np3grid),p1out8z( &0:np1,0:np2,0:np3),p2out8z(0:np1,0:np2,0:np3),p3out8z(0:np1,0:np2, &0:np3),p1cross(100),p2cross(100),p3cross(100),tout8z(0:nsave),uout &8z(0:np1,0:np2,0:np3,4*neqn,0:nsave),uout(0:np1,0:np2,0:np3,4,neqn &,0:nsave),xres8z(nxp8z),yres8z(nyp8z),zres8z(nzp8z),ures8z(neqn,nx &p8z,nyp8z,nzp8z) equivalence (uout,uout8z) allocatable iwrk8z(:),rwrk8z(:) C dimension iwrk8z(iiwk8z),rwrk8z(irwk8z) character*40 title logical linear,crankn,noupdt,nodist,fillin,evcmpx,adapt,plot,lsqfi &t,fdiff,solid,econ8z,ncon8z,restrt,gridid common/dtdp14/ sint8z(20),bint8z(20),slim8z(20),blim8z(20) common/dtdp15/ evlr8z,ev0r,evli8z,ev0i,evcmpx common/dtdp16/ p8z,evr8z(50),evi8z(50) common/dtdp19/ toler(neqnmx),adapt common/dtdp30/ econ8z,ncon8z common/dtdp45/ perdc(neqnmx) common/dtdp46/ eps8z,cgtl8z,npmx8z,itype common/dtdp52/ nxa8z,nya8z,nza8z,kd8z common/dtdp53/ work8z(nxp8z*nyp8z*nzp8z+9) common/dtdp64/ amin8z(4*neqnmx),amax8z(4*neqnmx) pi = 4.0*atan(1.d0) nxa8z = nxp8z nya8z = nyp8z nza8z = nzp8z kd8z = kdeg8z C############################################################################## C If you don't want to read the FINE PRINT, default NPROB. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you want to solve several similar problems in the same run, set +# C + NPROB equal to the number of problems you want to solve. Then NPROB +# C + loops through the main program will be done, with IPROB=1,...,NPROB, +# C + and you can make the problem parameters vary with IPROB. NPROB +# C + defaults to 1. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NPROB = 1 do 78755 iprob=1,nprob C############################################################################## C You may now define global parameters, which may be referenced in any # C of the "FORTRAN expressions" you input throughout the rest of this # C interactive session. You will be prompted alternately for parameter # C names and their values; enter a blank name when you are finished. # C # C Parameter names are valid FORTRAN variable names, starting in # C column 1. Thus each name consists of 1 to 6 alphanumeric characters, # C the first of which must be a letter. If the first letter is in the # C range I-N, the parameter must be an integer. # C # C Parameter values are either FORTRAN constants or FORTRAN expressions # C involving only constants and global parameters defined on earlier # C lines. They may also be functions of the problem number IPROB, if # C you are solving several similar problems in one run (NPROB > 1). Note # C that you are defining global CONSTANTS, not functions; e.g., parameter # C values may not reference any of the independent or dependent variables # C of your problem. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you define other parameters here later, using an editor, you must +# C + add them to COMMON block /PARM8Z/ everywhere this block appears, if +# C + they are to be "global" parameters. +# C + +# C + The variable PI is already included as a global parameter, with an +# C + accurate value 3.14159... +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- D = & 5 C############################################################################## C A collocation finite element method is used, with tri-cubic Hermite # C basis functions on the elements (small boxes) defined by the grid # C points: # C P1GRID(1),...,P1GRID(NP1GRID) # C P2GRID(1),...,P2GRID(NP2GRID) # C P3GRID(1),...,P3GRID(NP3GRID) # C You will first be prompted for NP1GRID, the number of P1-grid points, # C then for P1GRID(1),...,P1GRID(NP1GRID). Any points defaulted will be # C uniformly spaced between the points you define; the first and last # C points cannot be defaulted. Then you will be prompted similarly # C for the number and values of the P2 and P3-grid points. The limits # C on the parameters are then: # C P1GRID(1) < P1 < P1GRID(NP1GRID) # C P2GRID(1) < P2 < P2GRID(NP2GRID) # C P3GRID(1) < P3 < P3GRID(NP3GRID) # C # C############################################################################## call dtdpwx(p1grid,np1grid,0) call dtdpwx(p2grid,np2grid,0) call dtdpwx(p3grid,np3grid,0) C P1GRID DEFINED C-----------------------------------------> INPUT FROM GUI <------------------- P1GRID(1) = & -1 P1GRID(NP1GRID/2) = 0 C-----------------------------------------> INPUT FROM GUI <------------------- P1GRID(NP1GRID) = & 1 C P2GRID DEFINED C-----------------------------------------> INPUT FROM GUI <------------------- P2GRID(1) = & 0 C-----------------------------------------> INPUT FROM GUI <------------------- P2GRID(NP2GRID) = & 1 C P3GRID DEFINED C-----------------------------------------> INPUT FROM GUI <------------------- P3GRID(1) = & 0 C-----------------------------------------> INPUT FROM GUI <------------------- P3GRID(NP3GRID) = & 2*pi call dtdpwx(p1grid,np1grid,1) call dtdpwx(p2grid,np2grid,1) call dtdpwx(p3grid,np3grid,1) C############################################################################## C If you don't want to read the FINE PRINT, enter ISOLVE = 1. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + The following linear system solvers are available: +# C + +# C + 1. Sparse direct method +# C + Harwell Library routine MA27 (used by permission) is +# C + used to solve the (positive definite) "normal" +# C + equations A**T*A*x = A**T*b. The normal equations, +# C + which are essentially the equations which would result +# C + if a least squares finite element method were used +# C + instead of a collocation method, are substantially +# C + more ill-conditioned than the original system Ax = b, +# C + so it may be important to use high precision if this +# C + option is chosen. +# C + 2. Frontal method +# C + This is an out-of-core band solver. If you want to +# C + override the default number of rows in the buffer (11),+# C + set a new value for NPMX8Z in the main program. +# C + 3. Jacobi conjugate gradient iterative method +# C + A preconditioned conjugate gradient iterative method +# C + is used to solve the (positive definite) normal +# C + equations. High precision is also important if this +# C + option is chosen. (This solver is MPI-enhanced, if +# C + MPI is available.) If you want to override the +# C + default convergence tolerance, set a new relative +# C + tolerance CGTL8Z in the main program. +# C + 4. Local solver (normal equations) +# C + 5. Local solver (original equations) +# C + Choose these options ONLY if alterative linear system +# C + solvers have been installed locally. See subroutines +# C + (D)TD3M, (D)TD3N in file (d)subs.f for instructions +# C + on how to add local solvers. +# C + 6. MPI-based parallel band solver +# C + This is a parallel solver which runs efficiently on +# C + multiple processor machines, under MPI. It is a +# C + band solver, with the matrix distributed over the +# C + available processors. Choose this option ONLY if the +# C + solver has been activated locally. See subroutine +# C + (D)TD3O in file (d)subs.f for instructions on how to +# C + activate this solver and the MPI-enhancements to the +# C + conjugate gradient solver. +# C + +# C + Enter ISOLVE = 1,2,3,4,5 or 6 to select a linear system solver. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## ISOLVE = 1 C *******STEADY-STATE PROBLEM itype = 1 t0 = 0.0 dt = 1.0 crankn = .false. noupdt = .false. C############################################################################## C Is this a linear problem? ("linear" means all differential equations # C and all boundary conditions are linear) # C############################################################################## LINEAR = .TRUE. C Number of Newton iterations NSTEPS = 1 FDIFF = .FALSE. C############################################################################## C You may calculate one or more integrals (over the entire region) of # C some functions of the solution and its derivatives. How many integrals # C (NINT), if any, do you want to calculate? # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + In the FORTRAN program created by the preprocessor, the computed +# C + values of the integrals will be returned in the vector SINT8Z. If +# C + several iterations or time steps are done, only the last computed +# C + values are saved in SINT8Z (all values are printed). +# C + +# C + A limiting value, SLIM8Z(I), for the I-th integral can be set +# C + below in the main program. The computations will then stop +# C + gracefully whenever SINT8Z(I) > SLIM8Z(I), for any I=1...NINT. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- NINT = 1 C############################################################################## C You may calculate one or more boundary integrals (over the entire # C boundary) of some functions of the solution and its derivatives. How # C many boundary integrals (NBINT), if any, do you want to calculate? # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + In the FORTRAN program created by the preprocessor, the computed +# C + values of the integrals will be returned in the vector BINT8Z. If +# C + several iterations or time steps are done, only the last computed +# C + values are saved in BINT8Z (all values are printed). +# C + +# C + A limiting value, BLIM8Z(I), for the I-th boundary integral can be +# C + set below in the main program. The computations will then stop +# C + gracefully whenever BINT8Z(I) > BLIM8Z(I), for any I=1...NBINT. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NBINT = 1 lsqfit = .false. RESTRT = .FALSE. GRIDID = .TRUE. C############################################################################## C If you do not have any periodic boundary conditions, enter IPERDC=0. # C # C Enter IPERDC=1 for periodic conditions at P1 = P1GRID(1),P1GRID(NP1GRID)# C IPERDC=2 for periodic conditions at P2 = P2GRID(1),P2GRID(NP2GRID)# C IPERDC=3 for periodic conditions at P3 = P3GRID(1),P3GRID(NP3GRID)# C IPERDC=4 for periodic conditions on both P1 and P2 # C IPERDC=5 for periodic conditions on both P1 and P3 # C IPERDC=6 for periodic conditions on both P2 and P3 # C IPERDC=7 for periodic conditions on P1, P2 and P3. # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + When periodic boundary conditions are selected, they apply to all +# C + variables by default. To turn off periodic boundary conditions on +# C + the I-th variable, set PERDC(I) to 0 (or another appropriate value +# C + of IPERDC) below in the main program and set the desired boundary +# C + conditions in subroutine GB8Z, "by hand". +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- IPERDC = 3 C############################################################################## C The solution is saved on an NP1+1 by NP2+1 by NP3+1 rectangular grid # C covering the box (P1A,P1B) x (P2A,P2B) x (P3A,P3B). Enter values for # C P1A,P1B,P2A,P2B,P3A,P3B. These variables are usually defaulted. # C # C The defaults are P1A = P1GRID(1), P1B = P1GRID(NP1GRID) # C P2A = P2GRID(1), P2B = P2GRID(NP2GRID) # C P3A = P3GRID(1), P3B = P3GRID(NP3GRID) # C # C############################################################################## C defaults for p1a,p1b,p2a,p2b,p3a,p3b p1a = p1grid(1) p1b = p1grid(np1grid) p2a = p2grid(1) p2b = p2grid(np2grid) p3a = p3grid(1) p3b = p3grid(np3grid) C DEFINE P1A,P1B,P2A,P2B,P3A,P3B HERE: call dtdpx3(np1,np2,np3,p1a,p1b,p2a,p2b,p3a,p3b,hp18z,hp28z,hp38z, &p1out8z,p2out8z,p3out8z,npts8z) C SOLUTION SAVED EVERY NOUT ITERATIONS NOUT = NSTEPS C *******allocate workspace call dtdpqx(np1grid,np2grid,np3grid,isolve,neqn,ii8z,ir8z,iperdc) if (iiwk8z.gt.1) ii8z = iiwk8z if (irwk8z.gt.1) ir8z = irwk8z allocate (iwrk8z(ii8z),rwrk8z(ir8z)) C *******DRAW GRID LINES? PLOT = .TRUE. C *******call pde solver call dtdp3x(p1grid, p2grid, p3grid, np1grid,np2grid, np3grid, neqn &, p1out8z, p2out8z, p3out8z, uout,tout8z, npts8z, t0, dt, nsteps, &nout, nsave, crankn, noupdt, itype, linear, isolve, rwrk8z, ir8z, &iwrk8z, ii8z, iperdc, plot, lsqfit, fdiff, nint, nbint, restrt, gr &idid) deallocate (iwrk8z,rwrk8z) C *******read from restart file to array ures8z C call dtdpr3(1,xres8z,nxp8z,yres8z,nyp8z,zres8z,nzp8z,ures8z,neqn) C *******write array ures8z back to restart file C call dtdpr3(2,xres8z,nxp8z,yres8z,nyp8z,zres8z,nzp8z,ures8z,neqn) C *******call user-written postprocessor call postpr(tout8z,nsave,p1out8z,p2out8z,p3out8z,np1,np2,np3,uout, &neqn,p1grid,p2grid,p3grid,np1grid,np2grid,np3grid) C *******CROSS-SECTION CONTOUR PLOTS C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means C (possibly as modified by UPRINT,..) # C 2 Cx # C 3 Cy # C 4 Cz # C############################################################################## IVAR = 1 C P3 = CONSTANT CROSS-SECTIONS ics8z = 3 C############################################################################## C Plots of the variable or vector as a function of P1 and P2 will be # C made, at the output grid P3-points closest to # C P3 = P3CROSS(1),...,P3CROSS(NP3VALS) # C # C Enter values for NP3VALS and P3CROSS(1),...,P3CROSS(NP3VALS). # C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- NP3VALS = 2 P3CROSS(1) = & p3a+0*(p3b-p3a) P3CROSS(2) = & p3a+1*(p3b-p3a) C############################################################################## C If your region is rectangular, enter ITPLOT=0, and you need not read # C the FINE PRINT. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Enter: +# C + ITPLOT = 0 if you want a rectangular cross-section to be drawn, +# C + with axes P1 and P2 (recommended when ITRANS=0). +# C + ITPLOT = 1 if you want a polar cross-section to be drawn, with +# C + P1=radius, P2=angle (axes A1=P1*COS(P2) and +# C + A2=P1*SIN(P2); recommended when ITRANS=1,2,-2). +# C + ITPLOT = -1 if you want a polar cross-section to be drawn, with +# C + P2=radius, P1=angle (axes A1=P2*COS(P1) and +# C + A2=P2*SIN(P1); recommended when ITRANS=-1). +# C + ITPLOT = 2 if the desired cross-section is neither rectangular +# C + nor polar, and you want to define the axes A1,A2 +# C + yourself. (Sometimes recommended when ITRANS=3,-3.) +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## ITPLOT = 2 ical8z = 1 C############################################################################## C If you don't want to read the FINE PRINT, default A1MIN,A1MAX,A2MIN, # C and A2MAX. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + By default, at each cross-section the limits (A1MIN,A1MAX) and +# C + (A2MIN,A2MAX) on the axes A1 and A2 are chosen so that the plot just +# C + fills the plot space. For a non-rectangular region, this means the +# C + scale may vary from cross-section to cross-section. Enter values +# C + for A1MIN,A1MAX,A2MIN,A2MAX here if you want to ensure a common +# C + scale for all cross-sections. (You must specify all four variables, +# C + or default all four.) +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## A1MIN = 0.0 A1MAX = 0.0 A2MIN = 0.0 A2MAX = 0.0 ISET1 = 1 ISET2 = NSAVE ISINC = 1 C############################################################################## C If you don't want to read the FINE PRINT, enter 'no'. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Do you want to scale the axes on the plot so that the region is +# C + undistorted? Otherwise the axes will be scaled so that the figure +# C + approximately fills the plot space. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NODIST = .FALSE. C alow = amin8z(ivar) ahigh = amax8z(ivar) C############################################################################## C Enter lower (UMIN) and upper (UMAX) bounds for the contour values. UMIN # C and UMAX are often defaulted. # C # C Labeled contours will be drawn corresponding to the values # C # C UMIN + S*(UMAX-UMIN), for S=0.05,0.15,...0.95. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + By default, UMIN and UMAX are set to the minimum and maximum values +# C + of the variable to be plotted. For a common scaling, you may want +# C + to set UMIN=ALOW, UMAX=AHIGH. ALOW and AHIGH are the minimum and +# C + maximum values over all output points and over all saved time steps +# C + or iterations. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## UMIN = 0.0 UMAX = 0.0 C############################################################################## C Do you want two additional unlabeled contours to be drawn between each # C pair of labeled contours? # C############################################################################## FILLIN = .TRUE. C############################################################################## C Enter a title, WITHOUT quotation marks. A maximum of 40 characters # C are allowed. The default is no title. # C############################################################################## TITLE = ' ' TITLE = 'C ' call dtdprx(tout8z,nsave,iset1,iset2,isinc) do 78757 is8z=iset1,iset2,isinc do 78756 kzv8z=1,np3vals call dtdpzx(p3cross(kzv8z),p3a,p3b,np3,kz8z) call dtdpln(uout8z(0,0,0,ivar,is8z),np1,np2,np3,p1a,p1b,p2a,p2b,p3 &a,p3b,ics8z,ix8z,jy8z,kz8z,title,umin,umax,nodist,fillin,tout8z(is &8z),a1min,a1max,a2min,a2max,itplot,ical8z) 78756 continue 78757 continue C *******CROSS-SECTION SURFACE PLOTS C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means C (possibly as modified by UPRINT,..) # C 2 Cx # C 3 Cy # C 4 Cz # C############################################################################## IVAR = 1 C P3 = CONSTANT CROSS-SECTIONS ics8z = 3 C############################################################################## C Plots of the variable or vector as a function of P1 and P2 will be # C made, at the output grid P3-points closest to # C P3 = P3CROSS(1),...,P3CROSS(NP3VALS) # C # C Enter values for NP3VALS and P3CROSS(1),...,P3CROSS(NP3VALS). # C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- NP3VALS = 2 P3CROSS(1) = & p3a+0*(p3b-p3a) P3CROSS(2) = & p3a+1*(p3b-p3a) ISET1 = 1 ISET2 = NSAVE ISINC = 1 C############################################################################## C Enter the view latitude, VLAT, and the view longitude, VLON, desired # C for this plot, in degrees. VLAT and VLON must be between 10 and 80 # C degrees; each defaults to 45 degrees. VLAT and VLON are usually # C defaulted. # C############################################################################## VLON = 45.0 VLAT = 45.0 C alow = amin8z(ivar) ahigh = amax8z(ivar) C############################################################################## C Specify the range (UMIN,UMAX) for the dependent variable axis. UMIN # C and UMAX are often defaulted. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + By default, each plot will be scaled to just fit in the plot area. +# C + For a common scaling, you may want to set UMIN=ALOW, UMAX=AHIGH. +# C + ALOW and AHIGH are the minimum and maximum values over all output +# C + points and over all saved time steps or iterations. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## UMIN = 0.0 UMAX = 0.0 UMIN = & alow UMAX = & ahigh C############################################################################## C Enter a title, WITHOUT quotation marks. A maximum of 40 characters # C are allowed. The default is no title. # C############################################################################## TITLE = ' ' TITLE = 'C ' call dtdprx(tout8z,nsave,iset1,iset2,isinc) do 78759 is8z=iset1,iset2,isinc do 78758 kzv8z=1,np3vals call dtdpzx(p3cross(kzv8z),p3a,p3b,np3,kz8z) call dtdplo(p1out8z,p2out8z,p3out8z,uout8z(0,0,0,ivar,is8z),np1,np &2,np3,ics8z,ix8z,jy8z,kz8z,title,vlon,vlat,umin,umax,tout8z(is8z)) 78758 continue 78759 continue 78755 continue call endgks stop end subroutine tran8z(itrans,p1,p2,p3) implicit double precision (a-h,o-z) common /dtdp41/x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3,x11,x21,x31,x12,x2 &2,x32,x13,x23,x33,y11,y21,y31,y12,y22,y32,y13,y23,y33,z11,z21,z31, &z12,z22,z32,z13,z23,z33 common/parm8z/ pi,D C############################################################################## C You can solve problems in your region only if you can describe it by # C X = X(P1,P2,P3) # C Y = Y(P1,P2,P3) # C Z = Z(P1,P2,P3) # C with constant limits on the parameters P1,P2,P3. If your region is # C rectangular, enter ITRANS=0 and the trivial parameterization # C X = P1 # C Y = P2 # C Z = P3 # C will be used. Otherwise, you need to read the FINE PRINT below. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If P1,P2,P3 represent cylindrical, spherical or other non-Cartesian +# C + coordinates, you can reference the Cartesian coordinates X,Y,Z +# C + and derivatives of your unknowns with respect to these coordinates, +# C + when you define your PDE coefficients, boundary conditions, and +# C + volume and boundary integrals, if you enter ITRANS .NE. 0. Enter: +# C + ITRANS = 1, if P1,P2,P3 are cylindrical coordinates, that is, if +# C + P1=R, P2=Theta, P3=Z, where X = R*cos(Theta) +# C + Y = R*sin(Theta) +# C + Z = Z +# C + ITRANS = -1, same as ITRANS=1, but P1=Theta, P2=R, P3=Z +# C + ITRANS = 2, if P1,P2,P3 are spherical coordinates, that is, if +# C + P1=Rho, P2=Phi, P3=Theta, where +# C + X = Rho*sin(Phi)*cos(Theta) +# C + Y = Rho*sin(Phi)*sin(Theta) +# C + Z = Rho*cos(Phi) +# C + (Theta is longitude, Phi is measured from north pole) +# C + ITRANS = -2, same as ITRANS=2, but P1=Rho, P2=Theta, P3=Phi +# C + ITRANS = 3, to define your own coordinate transformation. In this +# C + case, you will be prompted to define X,Y,Z and their +# C + first and second derivatives in terms of P1,P2,P3. +# C + Because of symmetry, you will not be prompted for all +# C + of the second derivatives. If you make a mistake in +# C + computing any of these derivatives, PDE2D will usually +# C + be able to issue a warning message. (X1 = dX/dP1, etc) +# C + ITRANS = -3, same as ITRANS=3, but you will only be prompted to +# C + define X,Y,Z; their first and second derivatives will +# C + be approximated using finite differences. +# C + When ITRANS = -3 or 3, the first derivatives of X,Y,Z must all be +# C + continuous. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## ITRANS = -3 C-----------------------------------------> INPUT FROM GUI <------------------- X = & P1 C-----------------------------------------> INPUT FROM GUI <------------------- Y = & P2*R(P1)*cos(P3) C-----------------------------------------> INPUT FROM GUI <------------------- Z = & P2*R(P1)*sin(P3) return end subroutine pdes8z(yd8z,i8z,j8z,kint8z,p1,p2,p3,t,uu8z) implicit double precision (a-h,o-z) parameter (neqnmx= 99) C un8z(1,I),un8z(2,I),... hold the (rarely used) values C of UI,UI1,... from the previous iteration or time step common /dtdp5x/un8z(10,neqnmx) common /dtdp18/norm1,norm2,norm3 double precision norm1,norm2,norm3,normx,normy,normz dimension uu8z(10,neqnmx) common/parm8z/ pi,D C = uu8z(1, 1) C1 = uu8z(2, 1) C2 = uu8z(3, 1) C3 = uu8z(4, 1) C11= uu8z(5, 1) C22= uu8z(6, 1) C33= uu8z(7, 1) C12= uu8z(8, 1) C21= uu8z(8, 1) C13= uu8z(9, 1) C31= uu8z(9, 1) C23= uu8z(10, 1) C32= uu8z(10, 1) call dtdpcd(p1,p2,p3) call dtdpcb(p1,p2,p3,norm1,norm2,norm3,x,y,z,normx,normy,normz,3) call dtdpcc(p1,p2,p3, & C1,C2,C3,C11,C22,C33,C12,C13,C23, & x,y,z,Cx,Cy,Cz,Cxx,Cyy,Czz,Cxy,Cxz,Cyz, & Cyx,Czx,Czy,dvol,darea) Cnorm = Cx*normx + Cy*normy + Cz*normz if (i8z.eq.0) then yd8z = 0.0 C############################################################################## C Enter FORTRAN expressions for the functions whose integrals are to be # C calculated and printed. They may be functions of # C # C X,Y,Z,C,Cx,Cy,Cz,Cxx,Cyy,Czz,Cxy,Cxz,Cyz # C and (if applicable) T # C # C The parameters P1,P2,P3 and derivatives with respect to these may also # C be referenced (C1 = dC/dP1, etc): # C C1,C2,C3,C11,C22,C33,C12,C13,C23 # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you only want to integrate a function over part of the region, +# C + define that function to be zero in the rest of the region. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- C INTEGRAL DEFINED if (kint8z.eq. 1) yd8z = & C C############################################################################## C Enter FORTRAN expressions for the functions whose integrals are to be # C calculated and printed. They may be functions of # C # C X,Y,Z,C,Cx,Cy,Cz,Cxx,Cyy,Czz,Cxy,Cxz,Cyz # C and (if applicable) T # C # C The components (NORMx,NORMy,NORMz) of the unit outward normal vector # C may also be referenced. # C # C The parameters P1,P2,P3 and derivatives with respect to these may # C also be referenced: # C C1,C2,C3,C11,C22,C33,C12,C13,C23 # C You can also reference the normal derivative Cnorm. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you only want to integrate a function over part of the boundary, +# C + define that function to be zero on the rest of the boundary. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C BND. INTEGRAL1 DEFINED if (kint8z.eq.-1) yd8z = & D*Cnorm if (kint8z.gt.0) yd8z = yd8z*dvol if (kint8z.lt.0) yd8z = yd8z*darea else C############################################################################## C Now enter FORTRAN expressions to define the PDE coefficients, which # C may be functions of # C # C X,Y,Z,C,Cx,Cy,Cz,Cxx,Cyy,Czz,Cxy,Cxz,Cyz # C # C and, in some cases, of the parameter T. # C # C Recall that the PDE has the form # C # C F = 0 # C # C The parameters P1,P2,P3 and derivatives with respect to these may also # C be referenced (C1 = dC/dP1, etc): # C C1,C2,C3,C11,C22,C33,C12,C13,C23 # C############################################################################## if (j8z.eq.0) then yd8z = 0.0 C-----------------------------------------> INPUT FROM GUI <------------------- C F DEFINED if (i8z.eq. 1) yd8z = & D*(Cxx+Cyy+Czz) + Q(x) else endif endif return end function u8z(i8z,p1,p2,p3,t0) implicit double precision (a-h,o-z) common/parm8z/ pi,D call dtdpcd(p1,p2,p3) call dtdpcb(p1,p2,p3,zr18z,zr28z,zr38z,x,y,z,d18z,d28z,d38z,1) u8z = 0.0 return end subroutine gb8z(gd8z,ifac8z,i8z,j8z,p1,p2,p3,t,uu8z) implicit double precision (a-h,o-z) parameter (neqnmx= 99) dimension uu8z(10,neqnmx) C un8z(1,I),un8z(2,I),... hold the (rarely used) values C of UI,UI1,... from the previous iteration or time step common /dtdp5x/ un8z(10,neqnmx) common /dtdp18/norm1,norm2,norm3 double precision none,norm1,norm2,norm3,normx,normy,normz common/parm8z/ pi,D none = dtdplx(2) C = uu8z(1, 1) C1 = uu8z(2, 1) C2 = uu8z(3, 1) C3 = uu8z(4, 1) call dtdpcd(p1,p2,p3) call dtdpcb(p1,p2,p3,norm1,norm2,norm3,x,y,z,normx,normy,normz,3) call dtdpcb( & p1,p2,p3,C1,C2,C3,x,y,z,Cx,Cy,Cz,2) Cnorm = Cx*normx + Cy*normy + Cz*normz if (j8z.eq.0) gd8z = 0.0 C############################################################################## C Enter FORTRAN expressions to define the boundary condition functions, # C which may be functions of # C # C X,Y,Z,C,Cx,Cy,Cz and (if applicable) T # C # C Recall that the boundary conditions have the form # C # C G = 0 # C # C Enter NONE to indicate "no" boundary condition. # C # C The parameters P1,P2,P3 and derivatives with respect to these may also # C be referenced (C1 = dC/dP1, etc): # C C1,C2,C3 # C The components (NORMx,NORMy,NORMz) of the unit outward normal vector # C may also be referenced, as well as the normal derivative Cnorm. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If "no" boundary condition is specified, the PDE is enforced at +# C + points just inside the boundary (exactly on the boundary, if EPS8Z +# C + is set to 0 in the main program). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## if (ifac8z.eq. 1) then C############################################################################## C # C First define the boundary conditions on the face P1 = P1GRID(1). # C############################################################################## if (j8z.eq.0) then C-----------------------------------------> INPUT FROM GUI <------------------- C G DEFINED if (i8z.eq. 1) gd8z = & C else endif endif if (ifac8z.eq. 2) then C############################################################################## C # C Now define the boundary conditions on the face P1 = P1GRID(NP1GRID). # C############################################################################## if (j8z.eq.0) then C-----------------------------------------> INPUT FROM GUI <------------------- C G DEFINED if (i8z.eq. 1) gd8z = & none else endif endif if (ifac8z.eq. 3) then C############################################################################## C # C Now define the boundary conditions on the face P2 = P2GRID(1). # C############################################################################## if (j8z.eq.0) then C-----------------------------------------> INPUT FROM GUI <------------------- C G DEFINED if (i8z.eq. 1) gd8z = & none else endif endif if (ifac8z.eq. 4) then C############################################################################## C # C Now define the boundary conditions on the face P2 = P2GRID(NP2GRID). # C############################################################################## if (j8z.eq.0) then C-----------------------------------------> INPUT FROM GUI <------------------- C G DEFINED if (i8z.eq. 1) gd8z = & Cnorm else endif endif if (ifac8z.eq. 5) then C############################################################################## C # C Now define the boundary conditions on the face P3 = P3GRID(1). # C############################################################################## if (j8z.eq.0) then C PERIODIC BOUNDARY CONDITIONS SET (SEE IPERDC) C G DEFINED C if (i8z.eq. 1) gd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] else endif endif if (ifac8z.eq. 6) then C############################################################################## C # C Now define the boundary conditions on the face P3 = P3GRID(NP3GRID). # C############################################################################## if (j8z.eq.0) then C PERIODIC BOUNDARY CONDITIONS SET (SEE IPERDC) C G DEFINED C if (i8z.eq. 1) gd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] else endif endif return end subroutine pmod8z(p1,p2,p3,t,uu8z,uprint,uxprnt,uyprnt,uzprnt) implicit double precision (a-h,o-z) dimension uu8z(10,*),uprint(*),uxprnt(*),uyprnt(*),uzprnt(*) common/dtdp14/sint(20),bint(20),slim8z(20),blim8z(20) common/parm8z/ pi,D C = uu8z(1, 1) C1 = uu8z(2, 1) C2 = uu8z(3, 1) C3 = uu8z(4, 1) C11= uu8z(5, 1) C22= uu8z(6, 1) C33= uu8z(7, 1) C12= uu8z(8, 1) C21= uu8z(8, 1) C13= uu8z(9, 1) C31= uu8z(9, 1) C23= uu8z(10, 1) C32= uu8z(10, 1) call dtdpcd(p1,p2,p3) call dtdpcc(p1,p2,p3, & C1,C2,C3,C11,C22,C33,C12,C13,C23, & x,y,z,Cx,Cy,Cz,Cxx,Cyy,Czz,Cxy,Cxz,Cyz, & Cyx,Czx,Czy,dvol8z,dare8z) uxprnt( 1) = Cx uyprnt( 1) = Cy uzprnt( 1) = Cz C############################################################################## C If you don't want to read the FINE PRINT, default all of the following # C variables. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Normally, PDE2D saves the values of C,Cx,Cy,Cz at the output +# C + points. If different variables are to be saved (for later printing +# C + or plotting) the following functions can be used to re-define the +# C + output variables: +# C + define UPRINT(1) to replace C +# C + UXPRNT(1) Cx +# C + UYPRNT(1) Cy +# C + UZPRNT(1) Cz +# C + Each function may be a function of +# C + +# C + X,Y,Z,C,Cx,Cy,Cz,Cxx,Cyy,Czz,Cxy,Cxz,Cyz +# C + and (if applicable) T +# C + +# C + Each may also be a function of the integral estimates SINT(1),..., +# C + BINT(1),... +# C + +# C + The parameters P1,P2,P3 and derivatives with respect to these may +# C + also be referenced (C1 = dC/dP1, etc): +# C + C1,C2,C3,C11,C22,C33,C12,C13,C23 +# C + +# C + The default for each variable is no change, for example, UPRINT(1) +# C + defaults to C. Enter FORTRAN expressions for each of the +# C + following functions (or default). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C DEFINE UPRINT(*),UXPRNT(*),UYPRNT(*),UZPRNT(*) HERE: return end function axis8z(i8z,p1,p2,p3,ical8z) implicit double precision (a-h,o-z) common/parm8z/ pi,D call dtdpcd(p1,p2,p3) call dtdpcb(p1,p2,p3,zr18z,zr28z,zr38z,x,y,z,d18z,d28z,d38z,1) axis8z = 0.0 if (ical8z.eq. 1) then C############################################################################## C Define the axes A1 and A2 as functions of P1 and P2. A1 and A2 may # C also reference the constant P3 and the Cartesian coordinates X,Y,Z. # C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- C A1 DEFINED if (i8z.eq. 1) axis8z = & X C-----------------------------------------> INPUT FROM GUI <------------------- C A2 DEFINED if (i8z.eq. 2) axis8z = & sqrt(Y**2+Z**2) endif return end C dummy routines subroutine xy8z(i8z,iarc8z,s,x,y,s0,sf) implicit double precision (a-h,o-z) return end subroutine dis8z(x,y,ktri,triden,shape) implicit double precision (a-h,o-z) return end function fb8z(i8z,iarc8z,ktri,s,x,y,t) implicit double precision (a-h,o-z) fb8z = 0 return end subroutine postpr(tout,nsave,p1out,p2out,p3out,np1,np2,np3,uout,ne &qn,p1grid,p2grid,p3grid,np1grid,np2grid,np3grid) implicit double precision (a-h,o-z) dimension p1out(0:np1,0:np2,0:np3),p2out(0:np1,0:np2,0:np3) dimension p3out(0:np1,0:np2,0:np3),tout(0:nsave) dimension uout(0:np1,0:np2,0:np3,4,neqn,0:nsave) dimension p1grid(np1grid),p2grid(np2grid),p3grid(np3grid) common/parm8z/ pi,D common /dtdp27/ itask,npes,icomm common /dtdp46/ eps8z,cgtl8z,npmx8z,itype data lun,lud/0,47/ if (itask.gt.0) return C UOUT(I,J,K,IDER,IEQ,L) = U-sub-IEQ, if IDER=1 C Ux-sub-IEQ, if IDER=2 C Uy-sub-IEQ, if IDER=3 C Uz-sub-IEQ, if IDER=4 C (possibly as modified by UPRINT,..) C at the point (P1OUT(I,J,K) , P2OUT(I,J,K) , P3OUT(I,J,K)) C at time/iteration TOUT(L). C ******* ADD POSTPROCESSING CODE HERE: C IN THE EXAMPLE BELOW, MATLAB PLOTFILES pde2d.m, C pde2d.rdm CREATED (REMOVE COMMENTS TO ACTIVATE) if (lun.eq.0) then lun = 46 open (lun,file='pde2d.m') open (lud,file='pde2d.rdm') endif do 78750 l=0,nsave if (tout(l).ne.dtdplx(2)) nsave0 = l 78750 continue write (lud,78751) nsave0 write (lud,78751) neqn write (lud,78751) np1grid write (lud,78751) np2grid write (lud,78751) np3grid write (lud,78751) np1 write (lud,78751) np2 write (lud,78751) np3 78751 format (i8) do 78754 i=1,np1grid do 78753 j=1,np2grid do 78752 k=1,np3grid p1 = p1grid(i) p2 = p2grid(j) p3 = p3grid(k) call dtdpcd(p1,p2,p3) call dtdpcb(p1,p2,p3,z18z,z28z,z38z,x,y,z,d18z,d28z,d38z,1) write(lud,78764) x,y,z 78752 continue 78753 continue 78754 continue do 78757 i=0,np1 do 78756 j=0,np2 do 78755 k=0,np3 p1 = p1out(i,j,k) p2 = p2out(i,j,k) p3 = p3out(i,j,k) call dtdpcd(p1,p2,p3) call dtdpcb(p1,p2,p3,z18z,z28z,z38z,x,y,z,d18z,d28z,d38z,1) write(lud,78764) p1,p2,p3,x,y,z 78755 continue 78756 continue 78757 continue do 78763 l=0,nsave0 write (lud,78764) tout(l) do 78762 ieq=1,neqn do 78761 ider=1,4 do 78760 i=0,np1 do 78759 j=0,np2 do 78758 k=0,np3 write (lud,78764) uout(i,j,k,ider,ieq,l) 78758 continue 78759 continue 78760 continue 78761 continue 78762 continue 78763 continue 78764 format (e16.8) write (lun,*) '% Read solution from pde2d.rdm' write (lun,*) 'fid = fopen(''pde2d.rdm'');' write (lun,*) 'NSAVE = fscanf(fid,''%g'',1);' write (lun,*) 'NEQN = fscanf(fid,''%g'',1);' write (lun,*) 'NP1GRID = fscanf(fid,''%g'',1);' write (lun,*) 'NP2GRID = fscanf(fid,''%g'',1);' write (lun,*) 'NP3GRID = fscanf(fid,''%g'',1);' write (lun,*) 'NP1 = fscanf(fid,''%g'',1);' write (lun,*) 'NP2 = fscanf(fid,''%g'',1);' write (lun,*) 'NP3 = fscanf(fid,''%g'',1);' if (itype.eq.2) then write (lun,*) 'L0 = 0;' else write (lun,*) 'L0 = 1;' endif write (lun,*) 'T = zeros(NSAVE+1,1);' write (lun,*) 'P1 = zeros(NP2+1,NP1+1,NP3+1);' write (lun,*) 'P2 = zeros(NP2+1,NP1+1,NP3+1);' write (lun,*) 'P3 = zeros(NP2+1,NP1+1,NP3+1);' write (lun,*) 'U = zeros(NP2+1,NP1+1,NP3+1,NSAVE+1,4,NEQN);' write (lun,*) 'for JOB=1:2' write (lun,*) 'if (JOB==1)' write (lun,*) ' NX = NP1GRID;' write (lun,*) ' NY = NP2GRID;' write (lun,*) ' NZ = NP3GRID;' write (lun,*) 'else' write (lun,*) ' figure' write (lun,*) ' NX = NP1+1;' write (lun,*) ' NY = NP2+1;' write (lun,*) ' NZ = NP3+1;' write (lun,*) 'end' write (lun,*) 'X = zeros(NY,NX,NZ);' write (lun,*) 'Y = zeros(NY,NX,NZ);' write (lun,*) 'Z = zeros(NY,NX,NZ);' write (lun,*) 'for i=1:NX' write (lun,*) 'for j=1:NY' write (lun,*) 'for k=1:NZ' write (lun,*) ' if (JOB==2)' write (lun,*)' P1(j,i,k) = fscanf(fid,''%g'',1);' write (lun,*)' P2(j,i,k) = fscanf(fid,''%g'',1);' write (lun,*)' P3(j,i,k) = fscanf(fid,''%g'',1);' write (lun,*) ' end' write (lun,*)' X(j,i,k) = fscanf(fid,''%g'',1);' write (lun,*)' Y(j,i,k) = fscanf(fid,''%g'',1);' write (lun,*)' Z(j,i,k) = fscanf(fid,''%g'',1);' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'xmin = min(min(min(X(:,:,:))));' write (lun,*) 'xmax = max(max(max(X(:,:,:))));' write (lun,*) 'ymin = min(min(min(Y(:,:,:))));' write (lun,*) 'ymax = max(max(max(Y(:,:,:))));' write (lun,*) 'zmin = min(min(min(Z(:,:,:))));' write (lun,*) 'zmax = max(max(max(Z(:,:,:))));' write (lun,*) 'Xp = zeros(NY,NZ);' write (lun,*) 'Yp = zeros(NY,NZ);' write (lun,*) 'Zp = zeros(NY,NZ);' write (lun,*) 'Xp(:,:) = X(:,1,:);' write (lun,*) 'Yp(:,:) = Y(:,1,:);' write (lun,*) 'Zp(:,:) = Z(:,1,:);' write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')' write (lun,*) 'hold on' write (lun,*) 'Xp(:,:) = X(:,NX,:);' write (lun,*) 'Yp(:,:) = Y(:,NX,:);' write (lun,*) 'Zp(:,:) = Z(:,NX,:);' write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')' write (lun,*) 'hold on' write (lun,*) 'Xp = zeros(NX,NZ);' write (lun,*) 'Yp = zeros(NX,NZ);' write (lun,*) 'Zp = zeros(NX,NZ);' write (lun,*) 'Xp(:,:) = X(1,:,:);' write (lun,*) 'Yp(:,:) = Y(1,:,:);' write (lun,*) 'Zp(:,:) = Z(1,:,:);' write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')' write (lun,*) 'hold on' write (lun,*) 'Xp(:,:) = X(NY,:,:);' write (lun,*) 'Yp(:,:) = Y(NY,:,:);' write (lun,*) 'Zp(:,:) = Z(NY,:,:);' write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')' write (lun,*) 'hold on' write (lun,*) 'Xp = zeros(NY,NX);' write (lun,*) 'Yp = zeros(NY,NX);' write (lun,*) 'Zp = zeros(NY,NX);' write (lun,*) 'Xp(:,:) = X(:,:,1);' write (lun,*) 'Yp(:,:) = Y(:,:,1);' write (lun,*) 'Zp(:,:) = Z(:,:,1);' write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')' write (lun,*) 'hold on' write (lun,*) 'Xp(:,:) = X(:,:,NZ);' write (lun,*) 'Yp(:,:) = Y(:,:,NZ);' write (lun,*) 'Zp(:,:) = Z(:,:,NZ);' write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')' write (lun,*) 'axis([xmin xmax ymin ymax zmin zmax])' write (lun,*) 'xlabel(''X'')' write (lun,*) 'ylabel(''Y'')' write (lun,*) 'zlabel(''Z'')' write (lun,*) '%if (JOB==1)' write (lun,*) '% title(''Finite element grid'')' write (lun,*) '%else' write (lun,*) '% title(''Output grid'')' write (lun,*) '%end' write (lun,*) 'view(-15.0,30.0)' write (lun,*) 'end' write (lun,*) 'for l=0:NSAVE' write (lun,*) 'T(l+1) = fscanf(fid,''%g'',1);' write (lun,*) 'for ieq=1:NEQN' write (lun,*) 'for ider=1:4' write (lun,*) 'for i=0:NP1' write (lun,*) 'for j=0:NP2' write (lun,*) 'for k=0:NP3' write (lun,*) &' U(j+1,i+1,k+1,l+1,ider,ieq) = fscanf(fid,''%g'',1);' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) '% ******* Choose variable for cross-section plots' write (lun,*) '% (see comments in POSTPR)' write (lun,*) 'IEQ = 1;' write (lun,*) 'IDER = 1;' write (lun,*) '% Set ICS1=1 for P1=constant cross-sections' write (lun,*) 'ICS1 = 1;' write (lun,*) '% ICS2=1 for P2=constant cross-sections' write (lun,*) 'ICS2 = 1;' write (lun,*) '% ICS3=1 for P3=constant cross-sections' write (lun,*) 'ICS3 = 1;' write (lun,*) '% Number of cross-sections (>1)' write (lun,*) 'NCROSS = 5;' write (lun,*) 'umin = ', & 'min(min(min(min(U(:,:,:,L0+1:NSAVE+1,IDER,IEQ)))));' write (lun,*) 'umax = ', & 'max(max(max(max(U(:,:,:,L0+1:NSAVE+1,IDER,IEQ)))));' write (lun,*) 'for L=L0:NSAVE' write (lun,*) 'if(ICS1==1)' write (lun,*) 'Xp = zeros(NP2+1,NP3+1);' write (lun,*) 'Yp = zeros(NP2+1,NP3+1);' write (lun,*) 'Zp = zeros(NP2+1,NP3+1);' write (lun,*) 'Up = zeros(NP2+1,NP3+1);' write (lun,*) 'for I=0:max(1,fix(NP1/(NCROSS-1))):NP1' write (lun,*) ' figure' write (lun,*) ' Xp(:,:) = X(:,I+1,:);' write (lun,*) ' Yp(:,:) = Y(:,I+1,:);' write (lun,*) ' Zp(:,:) = Z(:,I+1,:);' write (lun,*) ' Up(:,:) = U(:,I+1,:,L+1,IDER,IEQ);' write (lun,*) ' surf(Xp,Yp,Zp,Up,''FaceColor'',''interp'')' write (lun,*) ' colorbar' write (lun,*) ' axis([xmin xmax ymin ymax zmin zmax])' write (lun,*) ' caxis([umin umax])' write (lun,*) ' xlabel(''X'')' write (lun,*) ' ylabel(''Y'')' write (lun,*) ' zlabel(''Z'')' write (lun,*) ' title(['' T = '',num2str(T(L+1)),'// & ''', P1 = '',num2str(P1(1,I+1,1))])' write (lun,*) ' view(-37.5,30.0)' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'if(ICS2==1)' write (lun,*) 'Xp = zeros(NP1+1,NP3+1);' write (lun,*) 'Yp = zeros(NP1+1,NP3+1);' write (lun,*) 'Zp = zeros(NP1+1,NP3+1);' write (lun,*) 'Up = zeros(NP1+1,NP3+1);' write (lun,*) 'for J=0:max(1,fix(NP2/(NCROSS-1))):NP2' write (lun,*) ' figure' write (lun,*) ' Xp(:,:) = X(J+1,:,:);' write (lun,*) ' Yp(:,:) = Y(J+1,:,:);' write (lun,*) ' Zp(:,:) = Z(J+1,:,:);' write (lun,*) ' Up(:,:) = U(J+1,:,:,L+1,IDER,IEQ);' write (lun,*) ' surf(Xp,Yp,Zp,Up,''FaceColor'',''interp'')' write (lun,*) ' colorbar' write (lun,*) ' axis([xmin xmax ymin ymax zmin zmax])' write (lun,*) ' caxis([umin umax])' write (lun,*) ' xlabel(''X'')' write (lun,*) ' ylabel(''Y'')' write (lun,*) ' zlabel(''Z'')' write (lun,*) ' title(['' T = '',num2str(T(L+1)),'// & ''', P2 = '',num2str(P2(J+1,1,1))])' write (lun,*) ' view(-37.5,30.0)' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'if(ICS3==1)' write (lun,*) 'Xp = zeros(NP2+1,NP1+1);' write (lun,*) 'Yp = zeros(NP2+1,NP1+1);' write (lun,*) 'Zp = zeros(NP2+1,NP1+1);' write (lun,*) 'Up = zeros(NP2+1,NP1+1);' write (lun,*) 'for K=0:max(1,fix(NP3/(NCROSS-1))):NP3' write (lun,*) ' figure' write (lun,*) ' Xp(:,:) = X(:,:,K+1);' write (lun,*) ' Yp(:,:) = Y(:,:,K+1);' write (lun,*) ' Zp(:,:) = Z(:,:,K+1);' write (lun,*) ' Up(:,:) = U(:,:,K+1,L+1,IDER,IEQ);' write (lun,*) ' surf(Xp,Yp,Zp,Up,''FaceColor'',''interp'')' write (lun,*) ' colorbar' write (lun,*) ' axis([xmin xmax ymin ymax zmin zmax])' write (lun,*) ' caxis([umin umax])' write (lun,*) ' xlabel(''X'')' write (lun,*) ' ylabel(''Y'')' write (lun,*) ' zlabel(''Z'')' write (lun,*) ' title(['' T = '',num2str(T(L+1)),'// & ''', P3 = '',num2str(P3(1,1,K+1))])' write (lun,*) ' view(-37.5,30.0)' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'p1mn = min(min(min(P1(:,:,:))));' write (lun,*) 'p1mx = max(max(max(P1(:,:,:))));' write (lun,*) 'p2mn = min(min(min(P2(:,:,:))));' write (lun,*) 'p2mx = max(max(max(P2(:,:,:))));' write (lun,*) 'p3mn = min(min(min(P3(:,:,:))));' write (lun,*) 'p3mx = max(max(max(P3(:,:,:))));' write (lun,*) 'for L=0:-1' write (lun,*) 'if(ICS1==1)' write (lun,*) 'Xp = zeros(NP2+1,NP3+1);' write (lun,*) 'Yp = zeros(NP2+1,NP3+1);' write (lun,*) 'Up = zeros(NP2+1,NP3+1);' write (lun,*) 'for I=0:max(1,fix(NP1/(NCROSS-1))):NP1' write (lun,*) ' figure' write (lun,*) ' Xp(:,:) = P2(:,I+1,:);' write (lun,*) ' Yp(:,:) = P3(:,I+1,:);' write (lun,*) ' Up(:,:) = U(:,I+1,:,L+1,IDER,IEQ);' write (lun,*) ' surf(Xp,Yp,Up)' write (lun,*) ' axis([p2mn p2mx p3mn p3mx umin umax])' write (lun,*) ' xlabel(''P2'')' write (lun,*) ' ylabel(''P3'')' write (lun,*) ' zlabel(''U'')' write (lun,*) ' title(['' T = '',num2str(T(L+1)),'// & ''', P1 = '',num2str(P1(1,I+1,1))])' write (lun,*) ' view(-37.5,30.0)' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'if(ICS2==1)' write (lun,*) 'Xp = zeros(NP1+1,NP3+1);' write (lun,*) 'Yp = zeros(NP1+1,NP3+1);' write (lun,*) 'Up = zeros(NP1+1,NP3+1);' write (lun,*) 'for J=0:max(1,fix(NP2/(NCROSS-1))):NP2' write (lun,*) ' figure' write (lun,*) ' Xp(:,:) = P1(J+1,:,:);' write (lun,*) ' Yp(:,:) = P3(J+1,:,:);' write (lun,*) ' Up(:,:) = U(J+1,:,:,L+1,IDER,IEQ);' write (lun,*) ' surf(Xp,Yp,Up)' write (lun,*) ' axis([p1mn p1mx p3mn p3mx umin umax])' write (lun,*) ' xlabel(''P1'')' write (lun,*) ' ylabel(''P3'')' write (lun,*) ' zlabel(''U'')' write (lun,*) ' title(['' T = '',num2str(T(L+1)),'// & ''', P2 = '',num2str(P2(J+1,1,1))])' write (lun,*) ' view(-37.5,30.0)' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'if(ICS3==1)' write (lun,*) 'Xp = zeros(NP2+1,NP1+1);' write (lun,*) 'Yp = zeros(NP2+1,NP1+1);' write (lun,*) 'Up = zeros(NP2+1,NP1+1);' write (lun,*) 'for K=0:max(1,fix(NP3/(NCROSS-1))):NP3' write (lun,*) ' figure' write (lun,*) ' Xp(:,:) = P1(:,:,K+1);' write (lun,*) ' Yp(:,:) = P2(:,:,K+1);' write (lun,*) ' Up(:,:) = U(:,:,K+1,L+1,IDER,IEQ);' write (lun,*) ' surf(Xp,Yp,Up)' write (lun,*) ' axis([p1mn p1mx p2mn p2mx umin umax])' write (lun,*) ' xlabel(''P1'')' write (lun,*) ' ylabel(''P2'')' write (lun,*) ' zlabel(''U'')' write (lun,*) ' title(['' T = '',num2str(T(L+1)),'// & ''', P3 = '',num2str(P3(1,1,K+1))])' write (lun,*) ' view(-37.5,30.0)' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' return end function R(x) implicit double precision (a-h,o-z) pi = 4.*atan(1.d0) if (x.le.0.0) then R = 1 else R = cos(pi/2*x) endif return end function Q(x) implicit double precision (a-h,o-z) if (x.le.0.0) then Q = 0 else Q = 1 endif return end