C **************************
C * PDE2D 9.1 MAIN PROGRAM *
C **************************
C##############################################################################
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ This is a steady-state elasticity problem: $#
C $ $#
C $ dS11/dX + dS12/dY + dS13/dZ = 0 $#
C $ dS12/dX + dS22/dY + dS23/dZ = 0 $#
C $ dS13/dX + dS23/dY + dS33/dZ = 0 $#
C $ where $#
C $ S11 = A*Ux + B*Vy + B*Wz S12 = C*(Uy+Vx) $#
C $ S22 = A*Vy + B*Ux + B*Wz S13 = C*(Uz+Wx) $#
C $ S33 = A*Wz + B*Ux + B*Vy S23 = C*(Vz+Wy) $#
C $ $#
C $ where (U,V,W) is the displacement vector and $#
C $ A = E*(1-Vnu)/(1+Vnu)/(1-2*Vnu) $#
C $ B = E*Vnu/(1+Vnu)/(1-2*Vnu) $#
C $ C = E/2.0/(1+Vnu) $#
C $ and E = 2, Vnu=0.33 are the elastic modulus and Poisson ratio, and $#
C $ Ux means derivative with respect to X. $#
C $ $#
C $ The object modeled is a torus, of major radius R0=5 and minor radius $#
C $ R1=4, with a smaller torus, of minor radius 0.2*R1, removed; ie, $#
C $ a cross-section of the torus is an annulus. $#
C $ $#
C $ Toroidal coordinates will be used: $#
C $ X = (R0 + P3*COS(P2))*COS(P1) $#
C $ Y = (R0 + P3*COS(P2))*SIN(P1) $#
C $ Z = P3*SIN(P2) $#
C $ where P1 is the (toroidal) angle measured within the major circle, $#
C $ P2 is the (poloidal) angle measured in the minor circle, and P3 is $#
C $ radial distance from the center of the minor circle. 0 < P1 < 2*PI, $#
C $ 0 < P2 < 2*PI, 0.2*R1 < P3 < R1. $#
C $ $#
C $ There are periodic boundary conditions on P1 and on P2, and at the $#
C $ inner surface (P3=0.2*R1), the displacements U,V and W are 0. On $#
C $ the outer surface (P3=R1), there is a unit inward boundary force, $#
C $ ie, the boundary force vector is -(NORMx,NORMy,NORMz), where $#
C $ (NORMx,NORMy,NORMz) is the unit outward normal to the boundary, in $#
C $ Cartesian coordinates. This means the boundary condition is: $#
C $ S11*NORMx + S12*NORMy + S13*NORMz = -NORMx $#
C $ S12*NORMx + S22*NORMy + S23*NORMz = -NORMy $#
C $ S13*NORMx + S23*NORMy + S33*NORMz = -NORMz $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: no $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
implicit real (a-h,o-z)
parameter (neqnmx= 99)
C##############################################################################
C NP1GRID = number of P1-grid lines #
C##############################################################################
PARAMETER (NP1GRID = 4)
C##############################################################################
C NP2GRID = number of P2-grid lines #
C##############################################################################
PARAMETER (NP2GRID = 8)
C##############################################################################
C NP3GRID = number of P3-grid lines #
C##############################################################################
PARAMETER (NP3GRID = 8)
C##############################################################################
C How many differential equations (NEQN) are there in your problem? #
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: NEQN = 3 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
PARAMETER (NEQN = 3)
C DIMENSIONS OF WORK ARRAYS
PARAMETER (IRWK8Z= 1241089)
PARAMETER (IIWK8Z= 2359297)
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: NP1 = 10 $#
C $ NP2 = 15 $#
C $ NP3 = 8 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
PARAMETER (NP1 = 10)
PARAMETER (NP2 = 15)
PARAMETER (NP3 = 8)
PARAMETER (NSAVE = 1)
common/parm8z/ pi,R0 ,R1 ,E ,VNU ,A ,B
&,C
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)
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/ tdp14/ sint8z(20),bint8z(20),slim8z(20),blim8z(20)
common/ tdp15/ evlr8z,ev0r,evli8z,ev0i,evcmpx
common/ tdp16/ p8z,evr8z(50),evi8z(50)
common/ tdp19/ toler(neqnmx),adapt
common/ tdp30/ econ8z,ncon8z
common/ tdp45/ perdc(neqnmx)
common/ tdp46/ eps8z,cgtl8z,npmx8z,itype
common/ tdp52/ nxa8z,nya8z,nza8z,kd8z
common/ tdp53/ work8z(nxp8z*nyp8z*nzp8z+9)
common/ tdp64/ 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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ press [RETURN] to default NPROB $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: R0 $#
C $ 5 $#
C $ R1 $#
C $ 4 $#
C $ E $#
C $ 2 $#
C $ VNU $#
C $ 0.33 $#
C $ A $#
C $ E*(1-VNU)/(1+VNU)/(1-2*VNU) $#
C $ B $#
C $ E*VNU/(1+VNU)/(1-2*VNU) $#
C $ C $#
C $ E/2.0/(1+VNU) $#
C $ [blank line] $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
R0 =
& 5
R1 =
& 4
E =
& 2
VNU =
& 0.33
A =
& E*(1-VNU)/(1+VNU)/(1-2*VNU)
B =
& E*VNU/(1+VNU)/(1-2*VNU)
C =
& E/2.0/(1+VNU)
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: NP1GRID = 4 $#
C $ P1GRID(1) = 0 $#
C $ P1GRID(NP1GRID) = 2*PI $#
C $ and default the other P1GRID points $#
C $ NP2GRID = 8 $#
C $ P2GRID(1) = 0 $#
C $ P2GRID(NP2GRID) = 2*PI $#
C $ and default the other P2GRID points $#
C $ NP3GRID = 8 $#
C $ P3GRID(1) = 0.2*R1 $#
C $ P3GRID(NP3GRID) = R1 $#
C $ and default the other P3GRID points $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
call tdpwx(p1grid,np1grid,0)
call tdpwx(p2grid,np2grid,0)
call tdpwx(p3grid,np3grid,0)
C P1GRID DEFINED
P1GRID(1) =
& 0
P1GRID(NP1GRID) =
& 2*PI
C P2GRID DEFINED
P2GRID(1) =
& 0
P2GRID(NP2GRID) =
& 2*PI
C P3GRID DEFINED
P3GRID(1) =
& 0.2*R1
P3GRID(NP3GRID) =
& R1
call tdpwx(p1grid,np1grid,1)
call tdpwx(p2grid,np2grid,1)
call tdpwx(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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: ISOLVE = 3 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
ISOLVE = 3
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: yes $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: NINT = 1 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: NBINT = 1 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: IPERDC = 4 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
IPERDC = 4
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ press [RETURN] to default P1A,P1B,P2A,P2B,P3A,P3B $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
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 tdpx3(np1,np2,np3,p1a,p1b,p2a,p2b,p3a,p3b,hp18z,hp28z,hp38z,
&p1out8z,p2out8z,p3out8z,npts8z)
C SOLUTION SAVED EVERY NOUT ITERATIONS
NOUT = NSTEPS
ii8z = iiwk8z
ir8z = irwk8z
C *******DRAW GRID LINES?
PLOT = .TRUE.
C *******call pde solver
call tdp3x(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)
C *******read from restart file to array ures8z
C call tdpr3(1,xres8z,nxp8z,yres8z,nyp8z,zres8z,nzp8z,ures8z,neqn)
C *******write array ures8z back to restart file
C call tdpr3(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)
C *******CROSS-SECTION VECTOR PLOTS
C P1 = CONSTANT CROSS-SECTIONS
ics8z = 1
C##############################################################################
C Plots of the variable or vector as a function of P2 and P3 will be #
C made, at the output grid P1-points closest to #
C P1 = P1CROSS(1),...,P1CROSS(NP1VALS) #
C #
C Enter values for NP1VALS and P1CROSS(1),...,P1CROSS(NP1VALS). #
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ We want a plot of the displacements at the cross-section P1 = 0, so $#
C $ enter: NP1VALS = 1 $#
C $ P1CROSS(1) = 0 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
NP1VALS = 1
P1CROSS(1) =
& 0
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 P2 and P3 (recommended when ITRANS=0,1,-1,2, +#
C + or -2) +#
C + ITPLOT = 1 if you want a polar cross-section to be drawn, with +#
C + P2=radius, P3=angle (axes A1=P2*COS(P3) and +#
C + A2=P2*SIN(P3)). +#
C + ITPLOT = -1 if you want a polar cross-section to be drawn, with +#
C + P3=radius, P2=angle (axes A1=P3*COS(P2) and +#
C + A2=P3*SIN(P2)). +#
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: ITPLOT = -1 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
ITPLOT = -1
ical8z = 1
A1MIN = 0.0
A1MAX = 0.0
A2MIN = 0.0
A2MAX = 0.0
C##############################################################################
C Enter values for IVAR1, IVAR2, IVAR3 to select the components Vr1,Vr2 #
C and Vr3 of the vector to be plotted. #
C IVAR1,IVAR2,IVAR3 = 1 means U (possibly as modified by UPRINT,...) #
C 2 Ux #
C 3 Uy #
C 4 Uz #
C 5 V #
C 6 Vx #
C 7 Vy #
C 8 Vz #
C 9 W #
C 10 Wx #
C 11 Wy #
C 12 Wz #
C #
C A vector plot of the in-plane components (Vr1,Vr2) will be made, with a #
C contour plot of the out-of-plane component Vr3 superimposed. #
C #
C Caution: If ITPLOT.NE.0, Vr1 and Vr2 are assumed to be the components #
C of the in-plane vector in the directions of the NEW axes A1,A2. #
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ At the cross-section p1=0, the in-plane components of the $#
C $ displacement field are U and W, and V is the out-of-plane component. $#
C $ Thus enter: $#
C $ enter: IVAR1 = 1 $#
C $ IVAR2 = 9 $#
C $ IVAR3 = 5 $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
IVAR1 = 1
IVAR2 = 9
IVAR3 = 5
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: no $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
NODIST = .FALSE.
C
a1mag = max(abs(amin8z(ivar1)),abs(amax8z(ivar1)))
a2mag = max(abs(amin8z(ivar2)),abs(amax8z(ivar2)))
C##############################################################################
C For the purpose of scaling the arrows, the ranges of the two components #
C of the vector are assumed to be (-VR1MAG,VR1MAG) and (-VR2MAG,VR2MAG). #
C Enter values for VR1MAG and VR2MAG. VR1MAG and VR2MAG are often #
C defaulted. #
C #
C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C + By default, VR1MAG and VR2MAG are the maxima of the absolute values +#
C + of the first and second components. For a common scaling, you may +#
C + want to set VR1MAG=A1MAG, VR2MAG=A2MAG. A1MAG, A2MAG are the +#
C + maxima of the absolute values over all output points and over all +#
C + saved time steps or iterations. +#
C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ press [RETURN] to default VR1MAG and VR2MAG $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
VR1MAG = 0.0
VR2MAG = 0.0
C
a3low = amin8z(ivar3)
a3high = amax8z(ivar3)
C##############################################################################
C Enter lower (VR3MIN) and upper (VR3MAX) bounds for the contour values #
C for the third component of the vector. VR3MIN and VR3MAX are often #
C defaulted. #
C #
C Contours will be drawn corresponding to the values #
C #
C VR3MIN + S*(VR3MAX-VR3MIN), for S=0.05,0.15,...0.95. #
C #
C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C + By default, VR3MIN and VR3MAX are set to the minimum and maximum +#
C + values of the variable to be plotted. For a common scaling, you may +#
C + want to set VR3MIN=A3LOW, VR3MAX=A3HIGH. A3LOW and A3HIGH are the +#
C + minimum and maximum values over all output points and over all saved +#
C + time steps or iterations. +#
C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ Press [RETURN] to default VR3MIN and VR3MAX $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
VR3MIN = 0.0
VR3MAX = 0.0
C##############################################################################
C Enter a title, WITHOUT quotation marks. A maximum of 40 characters #
C are allowed. The default is no title. #
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: Displacements at torus cross-section $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
TITLE = ' '
TITLE = 'Displacements at torus cross-section '
call tdprx(tout8z,nsave,iset1,iset2,isinc)
do 78757 is8z=iset1,iset2,isinc
do 78756 ixv8z=1,np1vals
call tdpzx(p1cross(ixv8z),p1a,p1b,np1,ix8z)
call tdplq(uout8z(0,0,0,ivar1,is8z),uout8z(0,0,0,ivar2,is8z),uout
&8z(0,0,0,ivar3,is8z),np1,np2,np3,p1a,p1b,p2a,p2b,p3a,p3b,ics8z,ix8
&z,jy8z,kz8z,title,vr1mag,vr2mag,vr3min,vr3max,nodist,tout8z(is8z),
&a1min,a1max,a2min,a2max,itplot,ical8z)
78756 continue
78757 continue
78755 continue
call endgks
stop
end
subroutine tran8z(itrans,p1,p2,p3)
implicit real (a-h,o-z)
common / tdp41/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,R0 ,R1 ,E ,VNU ,A ,B
&,C
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: ITRANS = -3 $#
C $ $#
C $ Then enter the following, when prompted: $#
C $ X = (R0 + P3*COS(P2))*COS(P1) $#
C $ Y = (R0 + P3*COS(P2))*SIN(P1) $#
C $ Z = P3*SIN(P2) $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
ITRANS = -3
X =
& (R0 + P3*COS(P2))*COS(P1)
Y =
& (R0 + P3*COS(P2))*SIN(P1)
Z =
& P3*SIN(P2)
return
end
subroutine pdes8z(yd8z,i8z,j8z,kint8z,p1,p2,p3,t,uu8z)
implicit real (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 / tdp5x/un8z(10,neqnmx)
common / tdp18/norm1,norm2,norm3
real norm1,norm2,norm3,normx,normy,normz
dimension uu8z(10,neqnmx)
common/parm8z/ pi,R0 ,R1 ,E ,VNU ,A ,B
&,C
U = uu8z(1, 1)
U1 = uu8z(2, 1)
U2 = uu8z(3, 1)
U3 = uu8z(4, 1)
U11= uu8z(5, 1)
U22= uu8z(6, 1)
U33= uu8z(7, 1)
U12= uu8z(8, 1)
U21= uu8z(8, 1)
U13= uu8z(9, 1)
U31= uu8z(9, 1)
U23= uu8z(10, 1)
U32= uu8z(10, 1)
V = uu8z(1, 2)
V1 = uu8z(2, 2)
V2 = uu8z(3, 2)
V3 = uu8z(4, 2)
V11= uu8z(5, 2)
V22= uu8z(6, 2)
V33= uu8z(7, 2)
V12= uu8z(8, 2)
V21= uu8z(8, 2)
V13= uu8z(9, 2)
V31= uu8z(9, 2)
V23= uu8z(10, 2)
V32= uu8z(10, 2)
W = uu8z(1, 3)
W1 = uu8z(2, 3)
W2 = uu8z(3, 3)
W3 = uu8z(4, 3)
W11= uu8z(5, 3)
W22= uu8z(6, 3)
W33= uu8z(7, 3)
W12= uu8z(8, 3)
W21= uu8z(8, 3)
W13= uu8z(9, 3)
W31= uu8z(9, 3)
W23= uu8z(10, 3)
W32= uu8z(10, 3)
call tdpcd(p1,p2,p3)
call tdpcb(p1,p2,p3,norm1,norm2,norm3,x,y,z,normx,normy,normz,3)
call tdpcc(p1,p2,p3,
& U1,U2,U3,U11,U22,U33,U12,U13,U23,
& x,y,z,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz,
& Uyx,Uzx,Uzy,dvol,darea)
Unorm = Ux*normx + Uy*normy + Uz*normz
call tdpcc(p1,p2,p3,
& V1,V2,V3,V11,V22,V33,V12,V13,V23,
& x,y,z,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz,
& Vyx,Vzx,Vzy,dvol,darea)
Vnorm = Vx*normx + Vy*normy + Vz*normz
call tdpcc(p1,p2,p3,
& W1,W2,W3,W11,W22,W33,W12,W13,W23,
& x,y,z,Wx,Wy,Wz,Wxx,Wyy,Wzz,Wxy,Wxz,Wyz,
& Wyx,Wzx,Wzy,dvol,darea)
Wnorm = Wx*normx + Wy*normy + Wz*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,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz #
C V,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz #
C W,Wx,Wy,Wz,Wxx,Wyy,Wzz,Wxy,Wxz,Wyz #
C and (if applicable) T #
C #
C The parameters P1,P2,P3 and derivatives with respect to these may also #
C be referenced (U1 = dU/dP1, etc): #
C U1,U2,U3,U11,U22,U33,U12,U13,U23 #
C V1,V2,V3,V11,V22,V33,V12,V13,V23 #
C W1,W2,W3,W11,W22,W33,W12,W13,W23 #
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ We want to compute the total volume change, which is the integral $#
C $ over our region of Ux+Vy+Wz ("true" value = -697.200). Thus, $#
C $ enter: INTEGRAL = Ux+Vy+Wz $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
C INTEGRAL DEFINED
if (kint8z.eq. 1) yd8z =
& Ux+Vy+Wz
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,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz #
C V,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz #
C W,Wx,Wy,Wz,Wxx,Wyy,Wzz,Wxy,Wxz,Wyz #
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 U1,U2,U3,U11,U22,U33,U12,U13,U23 #
C V1,V2,V3,V11,V22,V33,V12,V13,V23 #
C W1,W2,W3,W11,W22,W33,W12,W13,W23 #
C You can also reference the normal derivatives Unorm,Vnorm,Wnorm. #
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 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ By the divergence theorem, the volume integral of Ux+Vy+Wz should $#
C $ equal the boundary integral of U*NORMx+V*NORMy+W*NORMz, so we want $#
C $ to compute this boundary integral to verify the equality. Note that $#
C $ boundary integrals are computed over all 6 boundary pieces, but in $#
C $ this case the integrals at p1=0 and p1=2*pi will cancel, and $#
C $ similarly at p2=0 and p2=2*pi. Thus enter: $#
C $ BND. INTEGRAL = U*NORMx+V*NORMy+W*NORMz $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
C BND. INTEGRAL DEFINED
if (kint8z.eq. -1) yd8z =
& U*NORMx+V*NORMy+W*NORMz
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,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz #
C V,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz #
C W,Wx,Wy,Wz,Wxx,Wyy,Wzz,Wxy,Wxz,Wyz #
C #
C and, in some cases, of the parameter T. #
C #
C Recall that the PDEs have the form #
C #
C F1 = 0 #
C F2 = 0 #
C F3 = 0 #
C #
C The parameters P1,P2,P3 and derivatives with respect to these may also #
C be referenced (U1 = dU/dP1, etc): #
C U1,U2,U3,U11,U22,U33,U12,U13,U23 #
C V1,V2,V3,V11,V22,V33,V12,V13,V23 #
C W1,W2,W3,W11,W22,W33,W12,W13,W23 #
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ When asked if you want to write a FORTRAN block, $#
C $ enter: no $#
C $ then enter the following, when prompted: $#
C $ F1 = A*Uxx+B*Vyx+B*Wzx + C*(Uyy+Vxy) + C*(Uzz+Wxz) $#
C $ F2 = C*(Uyx+Vxx) + A*Vyy+B*Uxy+B*Wzy + C*(Vzz+Wyz) $#
C $ F3 = C*(Uzx+Wxx) + C*(Vzy+Wyy) + A*Wzz+B*Uxz+B*Vyz $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
if (j8z.eq.0) then
yd8z = 0.0
C F1 DEFINED
if (i8z.eq. 1) yd8z =
& A*Uxx+B*Vyx+B*Wzx + C*(Uyy+Vxy) + C*(Uzz+Wxz)
C F2 DEFINED
if (i8z.eq. 2) yd8z =
& C*(Uyx+Vxx) + A*Vyy+B*Uxy+B*Wzy + C*(Vzz+Wyz)
C F3 DEFINED
if (i8z.eq. 3) yd8z =
& C*(Uzx+Wxx) + C*(Vzy+Wyy) + A*Wzz+B*Uxz+B*Vyz
else
endif
endif
return
end
function u8z(i8z,p1,p2,p3,t0)
implicit real (a-h,o-z)
common/parm8z/ pi,R0 ,R1 ,E ,VNU ,A ,B
&,C
call tdpcd(p1,p2,p3)
call tdpcb(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 real (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 / tdp5x/ un8z(10,neqnmx)
common / tdp18/norm1,norm2,norm3
real none,norm1,norm2,norm3,normx,normy,normz
common/parm8z/ pi,R0 ,R1 ,E ,VNU ,A ,B
&,C
none = tdplx(2)
U = uu8z(1, 1)
U1 = uu8z(2, 1)
U2 = uu8z(3, 1)
U3 = uu8z(4, 1)
V = uu8z(1, 2)
V1 = uu8z(2, 2)
V2 = uu8z(3, 2)
V3 = uu8z(4, 2)
W = uu8z(1, 3)
W1 = uu8z(2, 3)
W2 = uu8z(3, 3)
W3 = uu8z(4, 3)
call tdpcd(p1,p2,p3)
call tdpcb(p1,p2,p3,norm1,norm2,norm3,x,y,z,normx,normy,normz,3)
call tdpcb(
& p1,p2,p3,U1,U2,U3,x,y,z,Ux,Uy,Uz,2)
Unorm = Ux*normx + Uy*normy + Uz*normz
call tdpcb(
& p1,p2,p3,V1,V2,V3,x,y,z,Vx,Vy,Vz,2)
Vnorm = Vx*normx + Vy*normy + Vz*normz
call tdpcb(
& p1,p2,p3,W1,W2,W3,x,y,z,Wx,Wy,Wz,2)
Wnorm = Wx*normx + Wy*normy + Wz*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,U,Ux,Uy,Uz, #
C V,Vx,Vy,Vz, #
C W,Wx,Wy,Wz and (if applicable) T #
C #
C Recall that the boundary conditions have the form #
C #
C G1 = 0 #
C G2 = 0 #
C G3 = 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 (U1 = dU/dP1, etc): #
C U1,U2,U3 #
C V1,V2,V3 #
C W1,W2,W3 #
C The components (NORMx,NORMy,NORMz) of the unit outward normal vector #
C may also be referenced, as well as the normal derivatives Unorm, #
C Vnorm,Wnorm. #
C #
C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C + If "no" boundary condition is specified, the corresponding PDE is +#
C + enforced at points just inside the boundary (exactly on the +#
C + boundary, if EPS8Z 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 PERIODIC BOUNDARY CONDITIONS SET (SEE IPERDC)
C G1 DEFINED
C if (i8z.eq. 1) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C G2 DEFINED
C if (i8z.eq. 2) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C G3 DEFINED
C if (i8z.eq. 3) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
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 PERIODIC BOUNDARY CONDITIONS SET (SEE IPERDC)
C G1 DEFINED
C if (i8z.eq. 1) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C G2 DEFINED
C if (i8z.eq. 2) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C G3 DEFINED
C if (i8z.eq. 3) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
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 PERIODIC BOUNDARY CONDITIONS SET (SEE IPERDC)
C G1 DEFINED
C if (i8z.eq. 1) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C G2 DEFINED
C if (i8z.eq. 2) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C G3 DEFINED
C if (i8z.eq. 3) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
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 PERIODIC BOUNDARY CONDITIONS SET (SEE IPERDC)
C G1 DEFINED
C if (i8z.eq. 1) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C G2 DEFINED
C if (i8z.eq. 2) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C G3 DEFINED
C if (i8z.eq. 3) gd8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
else
endif
endif
if (ifac8z.eq. 5) then
C##############################################################################
C #
C Now define the boundary conditions on the face P3 = P3GRID(1). #
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: G1 = U $#
C $ G2 = V $#
C $ G3 = W $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
if (j8z.eq.0) then
C G1 DEFINED
if (i8z.eq. 1) gd8z =
& U
C G2 DEFINED
if (i8z.eq. 2) gd8z =
& V
C G3 DEFINED
if (i8z.eq. 3) gd8z =
& W
else
endif
endif
if (ifac8z.eq. 6) then
C##############################################################################
C #
C Now define the boundary conditions on the face P3 = P3GRID(NP3GRID). #
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ enter: $#
C $ G1 = (A*Ux+B*Vy+B*Wz)*NORMx+C*(Uy+Vx)*NORMy+C*(Uz+Wx)*NORMz + NORMx $#
C $ G2 = C*(Uy+Vx)*NORMx+(A*Vy+B*Ux+B*Wz)*NORMy+C*(Vz+Wy)*NORMz + NORMy $#
C $ G3 = C*(Uz+Wx)*NORMx+C*(Vz+Wy)*NORMy+(A*Wz+B*Ux+B*Vy)*NORMz + NORMz $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
if (j8z.eq.0) then
C G1 DEFINED
if (i8z.eq. 1) gd8z =
& (A*Ux+B*Vy+B*Wz)*NORMx+C*(Uy+Vx)*NORMy+C*(Uz+Wx)*NORMz + NORMx
C G2 DEFINED
if (i8z.eq. 2) gd8z =
& C*(Uy+Vx)*NORMx+(A*Vy+B*Ux+B*Wz)*NORMy+C*(Vz+Wy)*NORMz + NORMy
C G3 DEFINED
if (i8z.eq. 3) gd8z =
& C*(Uz+Wx)*NORMx+C*(Vz+Wy)*NORMy+(A*Wz+B*Ux+B*Vy)*NORMz + NORMz
else
endif
endif
return
end
subroutine pmod8z(p1,p2,p3,t,uu8z,uprint,uxprnt,uyprnt,uzprnt)
implicit real (a-h,o-z)
dimension uu8z(10,*),uprint(*),uxprnt(*),uyprnt(*),uzprnt(*)
common/ tdp14/sint(20),bint(20),slim8z(20),blim8z(20)
common/parm8z/ pi,R0 ,R1 ,E ,VNU ,A ,B
&,C
U = uu8z(1, 1)
U1 = uu8z(2, 1)
U2 = uu8z(3, 1)
U3 = uu8z(4, 1)
U11= uu8z(5, 1)
U22= uu8z(6, 1)
U33= uu8z(7, 1)
U12= uu8z(8, 1)
U21= uu8z(8, 1)
U13= uu8z(9, 1)
U31= uu8z(9, 1)
U23= uu8z(10, 1)
U32= uu8z(10, 1)
V = uu8z(1, 2)
V1 = uu8z(2, 2)
V2 = uu8z(3, 2)
V3 = uu8z(4, 2)
V11= uu8z(5, 2)
V22= uu8z(6, 2)
V33= uu8z(7, 2)
V12= uu8z(8, 2)
V21= uu8z(8, 2)
V13= uu8z(9, 2)
V31= uu8z(9, 2)
V23= uu8z(10, 2)
V32= uu8z(10, 2)
W = uu8z(1, 3)
W1 = uu8z(2, 3)
W2 = uu8z(3, 3)
W3 = uu8z(4, 3)
W11= uu8z(5, 3)
W22= uu8z(6, 3)
W33= uu8z(7, 3)
W12= uu8z(8, 3)
W21= uu8z(8, 3)
W13= uu8z(9, 3)
W31= uu8z(9, 3)
W23= uu8z(10, 3)
W32= uu8z(10, 3)
call tdpcd(p1,p2,p3)
call tdpcc(p1,p2,p3,
& U1,U2,U3,U11,U22,U33,U12,U13,U23,
& x,y,z,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz,
& Uyx,Uzx,Uzy,dvol8z,dare8z)
uxprnt( 1) = Ux
uyprnt( 1) = Uy
uzprnt( 1) = Uz
call tdpcc(p1,p2,p3,
& V1,V2,V3,V11,V22,V33,V12,V13,V23,
& x,y,z,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz,
& Vyx,Vzx,Vzy,dvol8z,dare8z)
uxprnt( 2) = Vx
uyprnt( 2) = Vy
uzprnt( 2) = Vz
call tdpcc(p1,p2,p3,
& W1,W2,W3,W11,W22,W33,W12,W13,W23,
& x,y,z,Wx,Wy,Wz,Wxx,Wyy,Wzz,Wxy,Wxz,Wyz,
& Wyx,Wzx,Wzy,dvol8z,dare8z)
uxprnt( 3) = Wx
uyprnt( 3) = Wy
uzprnt( 3) = Wz
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 U,Ux,Uy,Uz,V,Vx, +#
C + Vy,Vz,W,Wx,Wy,Wz at the output points. If different +#
C + variables are to be saved (for later printing or plotting) the +#
C + following functions can be used to re-define the output variables: +#
C + define UPRINT(1) to replace U +#
C + UXPRNT(1) Ux +#
C + UYPRNT(1) Uy +#
C + UZPRNT(1) Uz +#
C + UPRINT(2) V +#
C + UXPRNT(2) Vx +#
C + UYPRNT(2) Vy +#
C + UZPRNT(2) Vz +#
C + UPRINT(3) W +#
C + UXPRNT(3) Wx +#
C + UYPRNT(3) Wy +#
C + UZPRNT(3) Wz +#
C + Each function may be a function of +#
C + +#
C + X,Y,Z,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz +#
C + V,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz +#
C + W,Wx,Wy,Wz,Wxx,Wyy,Wzz,Wxy,Wxz,Wyz +#
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 (U1 = dU/dP1, etc): +#
C + U1,U2,U3,U11,U22,U33,U12,U13,U23 +#
C + V1,V2,V3,V11,V22,V33,V12,V13,V23 +#
C + W1,W2,W3,W11,W22,W33,W12,W13,W23 +#
C + +#
C + The default for each variable is no change, for example, UPRINT(1) +#
C + defaults to U. Enter FORTRAN expressions for each of the +#
C + following functions (or default). +#
C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C $ press [RETURN] to default all output modification variables $#
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 9 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#
C##############################################################################
C DEFINE UPRINT(*),UXPRNT(*),UYPRNT(*),UZPRNT(*) HERE:
return
end
function axis8z(i8z,p1,p2,p3,ical8z)
implicit real (a-h,o-z)
common/parm8z/ pi,R0 ,R1 ,E ,VNU ,A ,B
&,C
call tdpcd(p1,p2,p3)
call tdpcb(p1,p2,p3,zr18z,zr28z,zr38z,x,y,z,d18z,d28z,d38z,1)
axis8z = 0.0
return
end
C dummy routines
subroutine xy8z(i8z,iarc8z,s,x,y,s0,sf)
implicit real (a-h,o-z)
return
end
subroutine dis8z(x,y,ktri,triden,shape)
implicit real (a-h,o-z)
return
end
function fb8z(i8z,iarc8z,ktri,s,x,y,t)
implicit real (a-h,o-z)
fb8z = 0
return
end
subroutine postpr(tout,nsave,p1out,p2out,p3out,np1,np2,np3,uout,ne
&qn)
implicit real (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)
common/parm8z/ pi,R0 ,R1 ,E ,VNU ,A ,B
&,C
common / tdp27/ itask,npes,icomm
common / tdp46/ 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 C! COMMENTS TO ACTIVATE)
C! if (lun.eq.0) then
C! lun = 46
C! open (lun,file='pde2d.m')
C! open (lud,file='pde2d.rdm')
C! write (lun,*) 'fid = fopen(''pde2d.rdm'');'
C! endif
C! do 78753 l=0,nsave
C! if (tout(l).ne. tdplx(2)) nsave0 = l
C!78753 continue
C! write (lud,78754) nsave0
C! write (lud,78754) neqn
C! write (lud,78754) np1
C! write (lud,78754) np2
C! write (lud,78754) np3
C!78754 format (i8)
C! do 78757 i=0,np1
C! do 78756 j=0,np2
C! do 78755 k=0,np3
C! p1 = p1out(i,j,k)
C! p2 = p2out(i,j,k)
C! p3 = p3out(i,j,k)
C! call tdpcd(p1,p2,p3)
C! call tdpcb(p1,p2,p3,z18z,z28z,z38z,x,y,z,d18z,d28z,d38z,1)
C! write(lud,78764) p1,p2,p3,x,y,z
C!78755 continue
C!78756 continue
C!78757 continue
C! do 78763 l=0,nsave0
C! write (lud,78764) tout(l)
C! do 78762 ieq=1,neqn
C! do 78761 ider=1,4
C! do 78760 i=0,np1
C! do 78759 j=0,np2
C! do 78758 k=0,np3
C! write (lud,78764) uout(i,j,k,ider,ieq,l)
C!78758 continue
C!78759 continue
C!78760 continue
C!78761 continue
C!78762 continue
C!78763 continue
C!78764 format (e16.8)
C! write (lun,*) 'NSAVE = fscanf(fid,''%g'',1);'
C! write (lun,*) 'NEQN = fscanf(fid,''%g'',1);'
C! write (lun,*) 'NP1 = fscanf(fid,''%g'',1);'
C! write (lun,*) 'NP2 = fscanf(fid,''%g'',1);'
C! write (lun,*) 'NP3 = fscanf(fid,''%g'',1);'
C! if (itype.eq.2) then
C! write (lun,*) 'L0 = 0;'
C! else
C! write (lun,*) 'L0 = 1;'
C! endif
C! write (lun,*) 'T = zeros(NSAVE+1,1);'
C! write (lun,*) 'P1 = zeros(NP2+1,NP1+1,NP3+1);'
C! write (lun,*) 'P2 = zeros(NP2+1,NP1+1,NP3+1);'
C! write (lun,*) 'P3 = zeros(NP2+1,NP1+1,NP3+1);'
C! write (lun,*) 'X = zeros(NP2+1,NP1+1,NP3+1);'
C! write (lun,*) 'Y = zeros(NP2+1,NP1+1,NP3+1);'
C! write (lun,*) 'Z = zeros(NP2+1,NP1+1,NP3+1);'
C! write (lun,*) 'U = zeros(NP2+1,NP1+1,NP3+1,NSAVE+1,4,NEQN);'
C! write (lun,*) '% Read solution from pde2d.rdm'
C! write (lun,*) 'for i=0:NP1'
C! write (lun,*) 'for j=0:NP2'
C! write (lun,*) 'for k=0:NP3'
C! write (lun,*)' P1(j+1,i+1,k+1) = fscanf(fid,''%g'',1);'
C! write (lun,*)' P2(j+1,i+1,k+1) = fscanf(fid,''%g'',1);'
C! write (lun,*)' P3(j+1,i+1,k+1) = fscanf(fid,''%g'',1);'
C! write (lun,*)' X(j+1,i+1,k+1) = fscanf(fid,''%g'',1);'
C! write (lun,*)' Y(j+1,i+1,k+1) = fscanf(fid,''%g'',1);'
C! write (lun,*)' Z(j+1,i+1,k+1) = fscanf(fid,''%g'',1);'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'for l=0:NSAVE'
C! write (lun,*) 'T(l+1) = fscanf(fid,''%g'',1);'
C! write (lun,*) 'for ieq=1:NEQN'
C! write (lun,*) 'for ider=1:4'
C! write (lun,*) 'for i=0:NP1'
C! write (lun,*) 'for j=0:NP2'
C! write (lun,*) 'for k=0:NP3'
C! write (lun,*)
C! &' U(j+1,i+1,k+1,l+1,ider,ieq) = fscanf(fid,''%g'',1);'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'xmin = min(min(min(X(:,:,:))));'
C! write (lun,*) 'xmax = max(max(max(X(:,:,:))));'
C! write (lun,*) 'ymin = min(min(min(Y(:,:,:))));'
C! write (lun,*) 'ymax = max(max(max(Y(:,:,:))));'
C! write (lun,*) 'zmin = min(min(min(Z(:,:,:))));'
C! write (lun,*) 'zmax = max(max(max(Z(:,:,:))));'
C! write (lun,*) 'p1mn = min(min(min(P1(:,:,:))));'
C! write (lun,*) 'p1mx = max(max(max(P1(:,:,:))));'
C! write (lun,*) 'p2mn = min(min(min(P2(:,:,:))));'
C! write (lun,*) 'p2mx = max(max(max(P2(:,:,:))));'
C! write (lun,*) 'p3mn = min(min(min(P3(:,:,:))));'
C! write (lun,*) 'p3mx = max(max(max(P3(:,:,:))));'
C! write (lun,*) '% Draw 3D region'
C! write (lun,*) 'Xp = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'Yp = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'Zp = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'Xp(:,:) = X(:,1,:);'
C! write (lun,*) 'Yp(:,:) = Y(:,1,:);'
C! write (lun,*) 'Zp(:,:) = Z(:,1,:);'
C! write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')'
C! write (lun,*) 'hold on'
C! write (lun,*) 'Xp(:,:) = X(:,NP1+1,:);'
C! write (lun,*) 'Yp(:,:) = Y(:,NP1+1,:);'
C! write (lun,*) 'Zp(:,:) = Z(:,NP1+1,:);'
C! write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')'
C! write (lun,*) 'hold on'
C! write (lun,*) 'Xp = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'Yp = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'Zp = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'Xp(:,:) = X(1,:,:);'
C! write (lun,*) 'Yp(:,:) = Y(1,:,:);'
C! write (lun,*) 'Zp(:,:) = Z(1,:,:);'
C! write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')'
C! write (lun,*) 'hold on'
C! write (lun,*) 'Xp(:,:) = X(NP2+1,:,:);'
C! write (lun,*) 'Yp(:,:) = Y(NP2+1,:,:);'
C! write (lun,*) 'Zp(:,:) = Z(NP2+1,:,:);'
C! write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')'
C! write (lun,*) 'hold on'
C! write (lun,*) 'Xp = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'Yp = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'Zp = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'Xp(:,:) = X(:,:,1);'
C! write (lun,*) 'Yp(:,:) = Y(:,:,1);'
C! write (lun,*) 'Zp(:,:) = Z(:,:,1);'
C! write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')'
C! write (lun,*) 'hold on'
C! write (lun,*) 'Xp(:,:) = X(:,:,NP3+1);'
C! write (lun,*) 'Yp(:,:) = Y(:,:,NP3+1);'
C! write (lun,*) 'Zp(:,:) = Z(:,:,NP3+1);'
C! write (lun,*) 'surf(Xp,Yp,Zp,''FaceColor'',''g'')'
C! write (lun,*) 'axis([xmin xmax ymin ymax zmin zmax])'
C! write (lun,*) 'xlabel(''X'')'
C! write (lun,*) 'ylabel(''Y'')'
C! write (lun,*) 'zlabel(''Z'')'
C! write (lun,*) 'title(''Region and output grid'')'
C! write (lun,*) '% Choose variable for cross-section plots'
C! write (lun,*) '% (see comments in POSTPR)'
C! write (lun,*) 'IEQ = 1;'
C! write (lun,*) 'IDER = 1;'
C! write (lun,*) '% Set ICS1=1 for P1=constant cross-sections'
C! write (lun,*) 'ICS1 = 1;'
C! write (lun,*) '% ICS2=1 for P2=constant cross-sections'
C! write (lun,*) 'ICS2 = 0;'
C! write (lun,*) '% ICS3=1 for P3=constant cross-sections'
C! write (lun,*) 'ICS3 = 0;'
C! write (lun,*) '% Number of cross-sections (>1)'
C! write (lun,*) 'NCROSS = 5;'
C! write (lun,*) 'umin = ',
C! & 'min(min(min(min(U(:,:,:,L0+1:NSAVE+1,IDER,IEQ)))));'
C! write (lun,*) 'umax = ',
C! & 'max(max(max(max(U(:,:,:,L0+1:NSAVE+1,IDER,IEQ)))));'
C! write (lun,*) 'for L=L0:NSAVE'
C! write (lun,*) 'if(ICS1==1)'
C! write (lun,*) 'Xp = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'Yp = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'Zp = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'Up = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'for I=0:max(1,fix(NP1/(NCROSS-1))):NP1'
C! write (lun,*) ' figure'
C! write (lun,*) ' Xp(:,:) = X(:,I+1,:);'
C! write (lun,*) ' Yp(:,:) = Y(:,I+1,:);'
C! write (lun,*) ' Zp(:,:) = Z(:,I+1,:);'
C! write (lun,*) ' Up(:,:) = U(:,I+1,:,L+1,IDER,IEQ);'
C! write (lun,*) ' surf(Xp,Yp,Zp,Up,''FaceColor'',''interp'')'
C! write (lun,*) ' colorbar'
C! write (lun,*) ' axis([xmin xmax ymin ymax zmin zmax])'
C! write (lun,*) ' caxis([umin umax])'
C! write (lun,*) ' xlabel(''X'')'
C! write (lun,*) ' ylabel(''Y'')'
C! write (lun,*) ' zlabel(''Z'')'
C! write (lun,*) ' title(['' T = '',num2str(T(L+1)),'//
C! & ''', P1 = '',num2str(P1(1,I+1,1))])'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'if(ICS2==1)'
C! write (lun,*) 'Xp = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'Yp = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'Zp = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'Up = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'for J=0:max(1,fix(NP2/(NCROSS-1))):NP2'
C! write (lun,*) ' figure'
C! write (lun,*) ' Xp(:,:) = X(J+1,:,:);'
C! write (lun,*) ' Yp(:,:) = Y(J+1,:,:);'
C! write (lun,*) ' Zp(:,:) = Z(J+1,:,:);'
C! write (lun,*) ' Up(:,:) = U(J+1,:,:,L+1,IDER,IEQ);'
C! write (lun,*) ' surf(Xp,Yp,Zp,Up,''FaceColor'',''interp'')'
C! write (lun,*) ' colorbar'
C! write (lun,*) ' axis([xmin xmax ymin ymax zmin zmax])'
C! write (lun,*) ' caxis([umin umax])'
C! write (lun,*) ' xlabel(''X'')'
C! write (lun,*) ' ylabel(''Y'')'
C! write (lun,*) ' zlabel(''Z'')'
C! write (lun,*) ' title(['' T = '',num2str(T(L+1)),'//
C! & ''', P2 = '',num2str(P2(J+1,1,1))])'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'if(ICS3==1)'
C! write (lun,*) 'Xp = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'Yp = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'Zp = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'Up = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'for K=0:max(1,fix(NP3/(NCROSS-1))):NP3'
C! write (lun,*) ' figure'
C! write (lun,*) ' Xp(:,:) = X(:,:,K+1);'
C! write (lun,*) ' Yp(:,:) = Y(:,:,K+1);'
C! write (lun,*) ' Zp(:,:) = Z(:,:,K+1);'
C! write (lun,*) ' Up(:,:) = U(:,:,K+1,L+1,IDER,IEQ);'
C! write (lun,*) ' surf(Xp,Yp,Zp,Up,''FaceColor'',''interp'')'
C! write (lun,*) ' colorbar'
C! write (lun,*) ' axis([xmin xmax ymin ymax zmin zmax])'
C! write (lun,*) ' caxis([umin umax])'
C! write (lun,*) ' xlabel(''X'')'
C! write (lun,*) ' ylabel(''Y'')'
C! write (lun,*) ' zlabel(''Z'')'
C! write (lun,*) ' title(['' T = '',num2str(T(L+1)),'//
C! & ''', P3 = '',num2str(P3(1,1,K+1))])'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'for L=L0:NSAVE'
C! write (lun,*) 'if(ICS1==1)'
C! write (lun,*) 'Xp = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'Yp = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'Up = zeros(NP2+1,NP3+1);'
C! write (lun,*) 'for I=0:max(1,fix(NP1/(NCROSS-1))):NP1'
C! write (lun,*) ' figure'
C! write (lun,*) ' Xp(:,:) = P2(:,I+1,:);'
C! write (lun,*) ' Yp(:,:) = P3(:,I+1,:);'
C! write (lun,*) ' Up(:,:) = U(:,I+1,:,L+1,IDER,IEQ);'
C! write (lun,*) ' surf(Xp,Yp,Up)'
C! write (lun,*) ' axis([p2mn p2mx p3mn p3mx umin umax])'
C! write (lun,*) ' xlabel(''P2'')'
C! write (lun,*) ' ylabel(''P3'')'
C! write (lun,*) ' zlabel(''U'')'
C! write (lun,*) ' title(['' T = '',num2str(T(L+1)),'//
C! & ''', P1 = '',num2str(P1(1,I+1,1))])'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'if(ICS2==1)'
C! write (lun,*) 'Xp = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'Yp = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'Up = zeros(NP1+1,NP3+1);'
C! write (lun,*) 'for J=0:max(1,fix(NP2/(NCROSS-1))):NP2'
C! write (lun,*) ' figure'
C! write (lun,*) ' Xp(:,:) = P1(J+1,:,:);'
C! write (lun,*) ' Yp(:,:) = P3(J+1,:,:);'
C! write (lun,*) ' Up(:,:) = U(J+1,:,:,L+1,IDER,IEQ);'
C! write (lun,*) ' surf(Xp,Yp,Up)'
C! write (lun,*) ' axis([p1mn p1mx p3mn p3mx umin umax])'
C! write (lun,*) ' xlabel(''P1'')'
C! write (lun,*) ' ylabel(''P3'')'
C! write (lun,*) ' zlabel(''U'')'
C! write (lun,*) ' title(['' T = '',num2str(T(L+1)),'//
C! & ''', P2 = '',num2str(P2(J+1,1,1))])'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'if(ICS3==1)'
C! write (lun,*) 'Xp = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'Yp = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'Up = zeros(NP2+1,NP1+1);'
C! write (lun,*) 'for K=0:max(1,fix(NP3/(NCROSS-1))):NP3'
C! write (lun,*) ' figure'
C! write (lun,*) ' Xp(:,:) = P1(:,:,K+1);'
C! write (lun,*) ' Yp(:,:) = P2(:,:,K+1);'
C! write (lun,*) ' Up(:,:) = U(:,:,K+1,L+1,IDER,IEQ);'
C! write (lun,*) ' surf(Xp,Yp,Up)'
C! write (lun,*) ' axis([p1mn p1mx p2mn p2mx umin umax])'
C! write (lun,*) ' xlabel(''P1'')'
C! write (lun,*) ' ylabel(''P2'')'
C! write (lun,*) ' zlabel(''U'')'
C! write (lun,*) ' title(['' T = '',num2str(T(L+1)),'//
C! & ''', P3 = '',num2str(P3(1,1,K+1))])'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
C! write (lun,*) 'end'
return
end