## Interactive session for prepared example 14

```
*******************************************************
****  Welcome to the PDE2D 9.6 Interactive Driver  ****
*******************************************************

PDE2D can solve systems of (linear or nonlinear) steady-state,
time-dependent and eigenvalue partial differential equations
in 1D intervals, general 2D regions, and in a wide range of simple
3D regions.  Ordinary differential equation systems can also be
solved.

problem.  The answers you give will be used to construct a PDE2D
FORTRAN program, which can then be compiled and linked with the PDE2D
runtime routines to produce an executable program.  The FORTRAN
driver program created will be well-documented and highly readable
(most of the interactive messages are repeated in the comments),
so that minor modifications or corrections can be made directly to
the FORTRAN program, without the need to work through a new
interactive session.

You can alternatively create your PDE2D FORTRAN program using the
PDE2D graphical user interface (GUI) ("pde2d_gui [progname]").  It
[RETURN]
is extraordinarily easy to set up problems using the PDE2D GUI,
which handles 0D and 1D problems, and 2D and 3D problems in "a wide
range of simple regions".  However, the PDE2D GUI cannot handle
complex regions, so if you have a complex 2D region you must use
this Interactive Driver.

If this is your first time to use PDE2D, you may want to work
through an example problem before trying one of your own.  Do you
want to work through a prepared example?
|---- Enter yes or no
yes
Several prepared examples are available.  Enter:

1 - to see a simple problem: a simply-supported elastic plate
equation, with a unit load concentrated at the midpoint of
a square.
2 - to see a more complex problem: a non-linear, steady-state
PDE, solved in an annulus.  Dirichlet (U = ...) boundary
conditions are imposed on part of the boundary, and Neumann
(dU/dn = ...) conditions are imposed on the other part, in
this example.  The initial triangulation is generated
automatically, and adaptive grid refinement is illustrated.
3 - to see an eigenvalue problem.  The region has a curved
interface across which material properties vary abruptly, in
this example.
4 - to see the first part of a thermal stress problem.  In this
part, the temperature distribution in a V-notched block is
calculated by solving the time-dependent heat conduction
equation, using adaptive time step control.
5 - to see the second part of a thermal stress problem.  In
this part, the stresses induced in the V-notched block by
thermal expansion are calculated, using the temperature
distribution output by example 4.  You must run example 4
[RETURN]
and save the tabular output before you can run example 5.
Examples 4 and 5 illustrate communication between problems.
6 - to see a 1D time-dependent integro-differential equation for
a financial math application.  In this problem there is a
term involving an integral of the solution, which requires
that we use PDE2D's feature for interpolating the solution
at the last saved time step, for use in the integral term.
7 - to see a waveguide problem (an eigenvalue problem in which
the eigenvalue appears nonlinearly).  This example shows how
to handle boundary conditions of different types on the same
arc, and how to produce a plot of a computed integral vs time.
8 - to see the Navier-Stokes equations (penalty formulation)
solved for a fluid flowing around a bend.
9 - to see a 3D elasticity problem, solved in a torus.  This
example illustrates the use of user-defined coordinate
transformations to handle more general 3D regions.
10 - to see a time-dependent wave equation (reduced to a system
of two PDEs), solved in a 3D box.
11 - to see a 3D eigenvalue problem (the Schrodinger equation
in a hydrogen atom).  This example illustrates the use of
spherical coordinates and periodic boundary conditions.
12 - to see a 3D eigenvalue problem, solved in a composite region
[RETURN]
consisting of two cylinders of different material properties.
13 - to see the axisymmetric Navier-Stokes equations solved in a
non-rectangular channel, using the collocation FEM.
14 - to see a 1D saturated/unsaturated water flow problem.
15 - to see a 1D version of the Schrodinger eigenvalue equation
of example 11.
0 - (no example)

If you select one of the examples, the correct answer for each
interactive question will be supplied after the question.
|---- Enter an integer value in the range 0 to 15
14
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ The example problem is                                               \$
\$                                                                      \$
\$     d/dt{Theta(H)} = d/dx{K(H)*(Hx+1)}         0 < X < 20            \$
\$                                                0 < T < 10            \$
\$ with  Hx+1 = 0    at X=0  (no flow)                                  \$
\$       H = 20      at X=20                                            \$
\$                                                                      \$
\$ and initial condition                                                \$
\$                                                                      \$
\$       H(X,0) = -20-X                                                 \$
\$                                                                      \$
\$ Here Theta(H) represents the moisture content, given by:             \$
\$      Theta(H) = Tr + B*(Ts-Tr)/(B+|H|**beta)     for H < 0           \$
\$               = Ts                               for H > 0           \$
\$      where Tr = 0.075, Ts = 0.287, B = 1,611,000, beta = 3.96        \$
\$ K(H) represents the hydraulic conductivity, given by:                \$
\$          K(H) =  Ks*A/(A+|H|**gamma)             for H < 0           \$
\$               =  Ks                              for H > 0           \$
\$      where Ks = 0.00944, A = 1,175,000, gamma = 4.74                 \$
\$ and H is the water pressure head.                                    \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
[RETURN]
[RETURN]
In what follows, when you are told to enter a 'FORTRAN expression',
this means any valid FORTRAN expression of 65 characters or less.
In this expression, you may include references to FORTRAN function
subprograms.  You may define these functions line by line at the end
of the interactive session, when prompted, or add them later using an
editor.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you enter a "#" in the first column of any input line, this       +
+ instructs the interactive driver to read this and subsequent input   +
+ lines from the file "pde2d.in".  A "#" in the first column of an     +
+ input line in the file "pde2d.in" (or an end-of-file) instructs the  +
+ driver to switch back to interactive input.                          +
+                                                                      +
+ All lines input during an interactive session are echo printed to    +
+ a file "echo.out".  You may want to modify this file and rename it   +
+ "pde2d.in", and read some or all of your input from this file during +
+ your next interactive session.                                       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

PDE2D can solve problems with 0,1,2 or 3 space variables.  Enter the
[RETURN]

0 - to solve a time-dependent ordinary differential equation system,
or an algebraic or algebraic eigenvalue system
1 - to solve problems in 1D intervals
2 - to solve problems in general 2D regions
3 - to solve problems in a wide range of simple 3D regions

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 1                                                             \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 0 to 3
1
Which finite element method do you want to use:

1. Galerkin method
2. Collocation method

Enter 1 or 2 to select a finite element method.
+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you have a problem with discontinuous material parameters, you    +
+ should use the Galerkin method.  If your partial differential        +
+ equations and boundary conditions are difficult to put into the      +
+ "divergence" form required by the Galerkin method, or if you have    +
+ periodic boundary conditions, use the collocation method.  The       +
+ collocation method produces approximate solutions with continuous    +
+ first derivatives; the Galerkin solutions are continuous but their   +
+ first derivatives are not.                                           +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 1                                                             \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 1 to 2
1
Is double precision mode to be used?  Double precision is recommended.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If double precision mode is used, variables and functions assigned   +
+ names beginning with a letter in the range A-H or O-Z will be DOUBLE +
+ PRECISION, and you should use double precision constants and FORTRAN +
+ expressions throughout; otherwise such variables and functions will  +
+ be of type REAL.  In either case, variables and functions assigned   +
+ names beginning with I,J,K,L,M or N will be of INTEGER type.         +
+                                                                      +
+ It is possible to convert a single precision PDE2D program to double +
+ precision after it has been created, using an editor.  Just change   +
+ all occurrences of "real" to "double precision"                      +
+                    " tdp" to "dtdp"  (note leading blank)            +
+ Any user-written code or routines must be converted "by hand", of    +
+ course.  To convert from double to single, reverse the changes.      +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
If you don't want to read the FINE PRINT, default NPROB.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you want to solve several similar problems in the same run, set   +
+ NPROB equal to the number of problems you want to solve.  Then NPROB +
+ loops through the main program will be done, with IPROB=1,...,NPROB, +
+ and you can make the problem parameters vary with IPROB.  NPROB      +
+ defaults to 1.                                                       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default NPROB                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NPROB =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

PDE2D solves the time-dependent system (note: U,A,F,FB,GB,U0 may be
vectors, C,RHO may be matrices):

C(X,T,U,Ux)*d(U)/dT = d/dX* A(X,T,U,Ux) - F(X,T,U,Ux)

d/dX* A(X,U,Ux) = F(X,U,Ux)

or the linear and homogeneous eigenvalue system:

d/dX* A(X,U,Ux) = F(X,U,Ux) + lambda*RHO(X)*U

with 'fixed' boundary conditions:

U = FB(X,[T])

or 'free' boundary conditions:

A*nx = GB(X,[T],U,Ux)      (nx=-1 at left endpoint, nx=+1 at right)

at two X values.
[RETURN]

For time-dependent problems there are also initial conditions:

U = U0(X)   at T=T0

Here Ux represents the (vector) function dU/dX.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If your PDEs involve the solution at points other than X, the        +
+ function                                                             +
+             (D)OLDSOL1(IDER,IEQ,XX,KDEG)                             +
+ will interpolate (using interpolation of degree KDEG=1,2 or 3) to XX +
+ the function saved in UOUT(*,IDER,IEQ,ISET) on the last time step or +
+ iteration (ISET) for which it has been saved.  Thus, for example, if +
+ IDER=1, this will return the latest value of component IEQ of the    +
+ solution at XX, assuming this has not been modified using UPRINT...  +
+ If your equations involve integrals of the solution, for example,    +
+ you can use (D)OLDSOL1 to approximate these using the solution from  +
+ the last time step or iteration.                                     +
+                                                                      +
+ CAUTION: For a steady-state or eigenvalue problem, you must reset    +
+ NOUT=1 if you want to save the solution each iteration.              +
[RETURN]
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ A system of NEQN complex partial differential equations must be      +
+ written as a system of 2*NEQN real equations, by separating the      +
+ equations into their real and imaginary parts.  However, note that   +
+ the complex arithmetic abilities of FORTRAN can be used to simplify  +
+ this separation.  For example, the complex PDE:                      +
+       I*Uxx = 1/(1+U**10),  where U = UR + UI*I                      +
+ would be difficult to split up analytically, but using FORTRAN       +
+ expressions it is easy:                                              +
+   A1 = -UIx,  F1 =  REAL(1.0/(1.0+CMPLX(UR,UI)**10))                 +
+   A2 =  URx,  F2 = AIMAG(1.0/(1.0+CMPLX(UR,UI)**10))                 +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
You may now define global parameters, which may be referenced in any
of the "FORTRAN expressions" you input throughout the rest of this
interactive session.  You will be prompted alternately for parameter
names and their values; enter a blank name when you are finished.

Parameter names are valid FORTRAN variable names, starting in
column 1.  Thus each name consists of 1 to 6 alphanumeric characters,
the first of which must be a letter.  If the first letter is in the
range I-N, the parameter must be an integer.

[RETURN]
Parameter values are either FORTRAN constants or FORTRAN expressions
involving only constants and global parameters defined on earlier
lines.  They may also be functions of the problem number IPROB, if
you are solving several similar problems in one run (NPROB > 1).  Note
that you are defining global CONSTANTS, not functions; e.g., parameter
values may not reference any of the independent or dependent variables

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you define other parameters here later, using an editor, you must +
+ add them to COMMON block /PARM8Z/ everywhere this block appears, if  +
+ they are to be "global" parameters.                                  +
+                                                                      +
+ The variable PI is already included as a global parameter, with an   +
+ accurate value 3.14159...                                            +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to skip global parameter definitions                  \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
Parameter name = (type blank line to terminate)

If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you want to be given a chance to write a FORTRAN block before the +
+ definitions of many functions?  If you answer 'no', you will still   +
+ be given a chance to write code before the definition of the PDE     +
+ coefficients, but not other functions.  Of course, you can always    +
+ add code later directly to the resulting program, using an editor.   +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
A Galerkin finite element method is used, with piecewise polynomial
basis functions on the subintervals defined by the grid points:
XGRID(1),XGRID(2),...,XGRID(NXGRID)
You will first be prompted for NXGRID, the number of X-grid points,
then for XGRID(1),...,XGRID(NXGRID).  Any points defaulted will be
uniformly spaced between the points you define; the first and last
points cannot be defaulted.  The interval over which the PDE system
is to be solved is then:
XGRID(1) < X < XGRID(NXGRID)

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NXGRID = 51                                                   \$
\$        XGRID(1) = 0.0                                                \$
\$        XGRID(NXGRID) = 20.0                                          \$
\$ and default XGRID(2),...,XGRID(NXGRID-1).                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NXGRID =
|---- Enter an integer value in the range 3 to +INFINITY
51
XGRID(1) =
|----Enter constant or FORTRAN expression-----------------------|
0.0
XGRID(2) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(3) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(4) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(5) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(6) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(7) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(8) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(9) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(10) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(11) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(12) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(13) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(14) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(15) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(16) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(17) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(18) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(19) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(20) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(21) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(22) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(23) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(24) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(25) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(26) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(27) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(28) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(29) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(30) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(31) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(32) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(33) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(34) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(35) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(36) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(37) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(38) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(39) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(40) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(41) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(42) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(43) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(44) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(45) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(46) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(47) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(48) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(49) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(50) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XGRID(NXGRID) =
|----Enter constant or FORTRAN expression-----------------------|
20.0
Enter the polynomial degree (1,2,3 or 4) desired.  A suggested value is
IDEG = 3.
+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ The spatial discretization error is O(h**2), O(h**3), O(h**4) or     +
+ O(h**5) when polynomials of degree 1,2,3 or 4, respectively, are     +
+ used.                                                                +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: IDEG = 3                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
IDEG =
|---- Enter an integer value in the range 1 to 4
3
What type of PDE problem do you want to solve?

2. a time-dependent problem
3. a linear, homogeneous eigenvalue problem

Enter 1,2 or 3 to select a problem type.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 2                                                             \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 1 to 3
2
Enter the initial time value (T0) and the final time value (TF), for
this time-dependent problem.  T0 defaults to 0.

TF is not required to be greater than T0.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: T0 = 0                                                        \$
\$        TF = 10                                                       \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
T0 =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
0
TF =
|----Enter constant or FORTRAN expression-----------------------|
10
Is this a linear problem? ("linear" means all differential equations
and all boundary conditions are linear).  If you aren't sure, it is
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Do you want the time step to be chosen adaptively?  If you answer
'yes', you will then be prompted to enter a value for TOLER(1), the
local relative time discretization error tolerance.  The default is
TOLER(1)=0.01.  If you answer 'no', a user-specified constant time step
will be used.  We suggest that you answer 'yes' and default TOLER(1)
(although for certain linear problems, a constant time step may be much
more efficient).

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If a negative value is specified for TOLER(1), then ABS(TOLER(1)) is +
+ taken to be the "absolute" error tolerance.  If a system of PDEs is  +
+ solved, by default the error tolerance specified in TOLER(1) applies +
+ to all variables, but the error tolerance for the J-th variable can  +
+ be set individually by specifying a value for TOLER(J) using an      +
+ editor, after the end of the interactive session.                    +
+                                                                      +
+ Each time step, two steps of size dt/2 are taken, and that solution  +
+ is compared with the result when one step of size dt is taken.  If   +
+ the maximum difference between the two answers is less than the      +
+ tolerance (for each variable), the time step dt is accepted (and the +
+ next step dt is doubled, if the agreement is "too" good); otherwise  +
+ dt is halved and the process is repeated.  Note that forcing the     +
[RETURN]
+ local (one-step) error to be less than the tolerance does not        +
+ guarantee that the global (cumulative) error is less than that value.+
+ However, as the tolerance is decreased, the global error should      +
+ decrease correspondingly.                                            +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$  then press [RETURN] to default TOLER(1)                             \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
TOLER(1) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

The time stepsize will be chosen adaptively, between an upper limit
of DTMAX = (TF-T0)/NSTEPS and a lower limit of 0.0001*DTMAX.  Enter
a value for NSTEPS (the minimum number of steps).

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you later turn off adaptive time step control, the time stepsize  +
+ will be constant, DT = (TF-T0)/NSTEPS.                               +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ So that the maximum time step will be (10.0-0.0)/20 = 0.5,           \$
\$ enter: NSTEPS = 20                                                   \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NSTEPS =
|----Enter constant or FORTRAN expression-----------------------|
20
If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Is the Crank-Nicolson scheme to be used to discretize time?  If you  +
+ answer 'no', a backward Euler scheme will be used.                   +
+                                                                      +
+ If a user-specified constant time step is chosen, the second order   +
+ Crank Nicolson method is recommended only for problems with very     +
+ well-behaved solutions, and the first order backward Euler scheme    +
+ should be used for more difficult problems.  In particular, do not   +
+ use the Crank Nicolson method if the left hand side of any PDE is    +
+ zero, for example, if a mixed elliptic/parabolic problem is solved.  +
+                                                                      +
+ If adaptive time step control is chosen, however, an extrapolation   +
+ is done between the 1-step and 2-step answers which makes the Euler  +
+ method second order, and the Crank-Nicolson method strongly stable.  +
+ Thus in this case, both methods have second order accuracy, and both +
+ are strongly stable.                                                 +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
[RETURN]
|---- Enter yes or no
yes
How many differential equations (NEQN) are there in your problem?
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NEQN = 1                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NEQN =
|---- Enter an integer value in the range 1 to 99
1
You may now choose names for the component(s) of the (possibly vector)
solution U.  Each must be an alphanumeric string of one to three
characters, beginning with a letter in the range A-H or O-Z.  The
variable names X,T and A must not be used.  The name should start
in column 1.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: U1 = H                                                        \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
U1 =
H
If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Is this problem symmetric?  The problem is called symmetric if the   +
+ Jacobian matrices                                                    +
+                                                                      +
+          F1.H     F1.Hx  ...                                         +
+          A1.H     A1.Hx  ...                                         +
+            .          .                                              +
+            .          .                                              +
+ and (if NEQN > 1)                                                    +
+          C11        C12      ...                                     +
+          C21        C22      ...                                     +
+            .          .                                              +
+            .          .                                              +
+ are always symmetric, where F1.H means d(F1)/d(H), and similarly     +
+ for the other terms.  The Jacobian of GB must also be symmetric,     +
+ that is, d(GB_i)/d(U_j) = d(GB_j)/d(U_i), and GB_i must not depend   +
+ on Hx,...                                                            +
+                                                                      +
+ The memory and execution time are decreased if the problem is known  +
+ to be symmetric.  If you answer 'yes' and your problem is not        +
[RETURN]
+ symmetric, you will usually get a warning message, but you may get   +
+ incorrect answers with no warning, so only answer 'yes' if you are   +
+ sure.                                                                +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

If you don't want to read the FINE PRINT, enter 'no'.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
If you don't want to read the FINE PRINT, enter 'yes' (strongly
recommended).

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ The partial derivatives of some of the PDE and boundary condition    +
+ coefficients are required by PDE2D.  These may be calculated         +
+ automatically using a finite difference approximation, or supplied   +
+ by the user.  Do you want them to be calculated automatically?       +
+                                                                      +
+ If you answer 'yes', you will not be asked to supply the derivatives,+
+ but there is a small risk that the inaccuracies introduced by the    +
+ finite difference approximation may cause the Newton iteration       +
+ to converge more slowly or to diverge, especially if low precision   +
+ is used.  This risk is very low, however, and since answering 'no'   +
+ means you may have to compute many partial derivatives, it is        +
+ recommended you answer 'yes' unless you have some reason to believe  +
+ there is a problem with the finite difference approximations.        +
+                                                                      +
+ If you supply analytic partial derivatives, PDE2D will do some spot  +
+ checking and can usually issue a warning if any are supplied         +
+ incorrectly.                                                         +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
[RETURN]
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
You may calculate one or more integrals (over the entire region) of
some functions of the solution and its derivatives.  How many integrals
(NINT), if any, do you want to calculate?

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ In the FORTRAN program created by the preprocessor, the computed     +
+ values of the integrals will be returned in the vector SINT8Z.  If   +
+ several iterations or time steps are done, only the last computed    +
+ values are saved in SINT8Z (all values are printed).                 +
+                                                                      +
+ A limiting value, SLIM8Z(I), for the I-th integral can be set        +
+ below in the main program.  The computations will then stop          +
+ gracefully whenever SINT8Z(I) > SLIM8Z(I), for any I=1...NINT.       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NINT = 1                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NINT =
|---- Enter an integer value in the range 0 to 20
1
Enter FORTRAN expressions for the functions whose integrals are to be
calculated and printed.  They may be functions of

X,H,Hx and (if applicable) T

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you only want to integrate a function over part of the interval,  +
+ define that function to be zero on the rest of the interval.         +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ We want to calculate the total moisture content, so                  \$
\$ enter: INTEGRAL = THETA(H,0)                                         \$
\$ The FORTRAN function THETA(H,0) (=Theta(H)) will be defined later.   \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
INTEGRAL =
|----Enter constant or FORTRAN expression-----------------------|
THETA(H,0)
You may calculate one or more boundary integrals (over the entire
boundary) of some functions of the solution and its derivatives.  How
many boundary integrals (NBINT), if any, do you want to calculate?

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ In the FORTRAN program created by the preprocessor, the computed     +
+ values of the integrals will be returned in the vector BINT8Z.  If   +
+ several iterations or time steps are done, only the last computed    +
+ values are saved in BINT8Z (all values are printed).                 +
+                                                                      +
+ A limiting value, BLIM8Z(I), for the I-th boundary integral can be   +
+ set below in the main program.  The computations will then stop      +
+ gracefully whenever BINT8Z(I) > BLIM8Z(I), for any I=1...NBINT.      +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NBINT = 1                                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NBINT =
|---- Enter an integer value in the range 0 to 20
1
Enter FORTRAN expressions for the functions whose "integrals" (sum
over two boundary points) are to be calculated and printed.  They may
be functions of

X,H,Hx and (if applicable) T

The unit outward normal, NORMx (=1 at right end point, -1 at left),
may also be referenced.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you only want to "integrate" a function over one boundary point,  +
+ define that function to be zero at the other point.                  +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ We want to calculate the flux, K(H)*(Hx+1), across the boundary, so  \$
\$ enter: BND. INTEGRAL = RK(H)*(Hx+1)*NORMx                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
BND. INTEGRAL =
|----Enter constant or FORTRAN expression-----------------------|
RK(H)*(Hx+1)*NORMx
If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you have point source terms in your PDEs, involving "Dirac Delta" +
+ functions?  A Dirac Delta function DEL(x-xd0) is a function whose    +
+ integral is 1, but which is zero everywhere except at a singular     +
+ point xd0.                                                           +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Now enter FORTRAN expressions to define the PDE coefficients, which
may be functions of

X,T,H,Hx

Recall that the PDE has the form

C*d(H)/dT = d/dX*A - F

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ The partial differential equation may be written                     \$
\$                                                                      \$
\$     d(Theta)/dH*Ht = d/dx{K(H)*(Hx+1)}                               \$
\$                                                                      \$
\$ When asked if you want to write a FORTRAN block,                     \$
\$ enter: no                                                            \$
\$ then enter the following, when prompted:                             \$
\$    C = THETA(H,1)       F = 0          A = RK(H)*(Hx+1)              \$
\$ The FORTRAN functions THETA(H,1) (=d(Theta)/dH) and RK(H) (=K(H))    \$
\$ will be defined later.                                               \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$

Do you want to write a FORTRAN block to define some parameters to be
used in the definition of these coefficients?
|---- Enter yes or no
no
C =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
THETA(H,1)
F =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
0
A =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
RK(H)*(Hx+1)
Now the initial values must be defined using FORTRAN expressions.
They may be functions of X, and may also reference the initial time T0.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: H0 = -20 - X                                                  \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
H0 =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
-20 - X
If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you want to read the initial conditions from the restart file,    +
+ if it exists (and use the conditions supplied above if it does not   +
+ exist)?                                                              +
+                                                                      +
+ If so, PDE2D will dump the final solution at the end of each run     +
+ into a restart file "pde2d.res".  Thus the usual procedure for       +
+ using this dump/restart option is to make sure there is no restart   +
+ file in your directory left over from a previous job, then the       +
+ first time you run this job, the initial conditions supplied above   +
+ will be used, but on the second and subsequent runs the restart file +
+ from the previous run will be used to define the initial conditions. +
+                                                                      +
+ You can do all the "runs" in one program, by setting NPROB > 1.      +
+ Each pass through the DO loop, T0,TF,NSTEPS and possibly other       +
+ parameters may be varied, by making them functions of IPROB.         +
+                                                                      +
+ If the 2D or 3D collocation method is used, the coordinate           +
+ transformation should not change between dump and restart.           +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
[RETURN]
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Do you have 'fixed' boundary conditions at the LEFT endpoint, that
is, are you going to specify values for all unknowns there?
If you answer 'no', you can specify boundary conditions of the
more general 'free' type.
+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If some, but not all, of the unknowns are to be specified, you       +
+ should specify 'free' boundary conditions, since even 'fixed'        +
+ type conditions of the form:                                         +
+                   Ui = FBi(X,[T])                                    +
+ can be expressed as 'free' boundary conditions in the form:          +
+               Ai*nx = zero(Ui-FBi(X,[T]))                            +
+ where zero(f) = Big_Number*f is a PDE2D-supplied function.           +
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Note: while 'fixed' boundary conditions are enforced exactly, 'free' +
+ boundary conditions are not; the greater the overall solution        +
+ accuracy, the more closely they are satisfied.                       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Do you have 'fixed' boundary conditions at the RIGHT endpoint, that
is, are you going to specify values for all unknowns there?
If you answer 'no', you can specify boundary conditions of the
more general 'free' type.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes

Define boundary conditions at the RIGHT endpoint X = XGRID(NXGRID).
Enter a FORTRAN expression to define FB at this endpoint.  It may
be a function of X and (if applicable) T.

Recall that fixed boundary conditions have the form

H = FB

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ The right endpoint boundary condition is H = 20, so                  \$
\$ enter: FB = 20                                                       \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
FB =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
20

Define boundary conditions at the LEFT endpoint X = XGRID(1).
[RETURN]
Enter a FORTRAN expression to define GB at this endpoint.  It may be
a function of

X,H,Hx and (if applicable) T

Recall that free boundary conditions have the form

A*nx = GB       (nx=-1 at left endpoint, nx=+1 at right)

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If A=0 (eg, if the PDE is first order), then setting GB=0 is         +
+ equivalent to setting "no" boundary condition (which is sometimes    +
+ appropriate for first order PDEs), because the boundary condition    +
+ A*nx = GB reduces to 0=0.                                            +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ The left endpoint boundary condition is Hx+1 = 0, or A = 0, so       \$
\$ enter: GB = 0                                                        \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
GB =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
0
If you don't want to read the FINE PRINT, default all of the following
variables.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Normally, PDE2D saves the values of H,A at the output points.        +
+ If different variables are to be saved (for later printing or        +
+ plotting) the following functions can be used to re-define the       +
+ output variables:                                                    +
+    define UPRINT(1) to replace H                                     +
+           APRINT(1)            A                                     +
+ Each function may be a function of                                   +
+                                                                      +
+           X,H,Hx,A and (if applicable) T                             +
+                                                                      +
+ Each may also be a function of the initial triangle number KTRI and  +
+ the integral estimates SINT(1),...,BINT(1),...                       +
+                                                                      +
+ The default for each variable is no change, for example, UPRINT(1)   +
+ defaults to H.  Enter FORTRAN expressions for each of the            +
+ following functions (or default).                                    +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
[RETURN]
\$ press [RETURN] to default all output modification variables          \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
Replace H for postprocessing?
UPRINT(1) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

Replace A for postprocessing?
APRINT(1) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

The solution is saved on a uniform grid of NX+1 points
XA + I*(XB-XA)/NX
I=0,...,NX.  Enter a value for NX (suggested value = 50).

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you want to save the solution at an arbitrary user-specified set  +
+ of NX+1 points, enter -NX.                                           +
+                                                                      +
+ If you set NEAR8Z=1 in the main program, the values saved at each    +
+ output point will actually be the solution as evaluated at a nearby  +
+ collocation or integration point.  For most problems this obviously  +
+ will produce less accurate output or plots, but for certain (rare)   +
+ problems, a solution component may be much less noisy when plotted   +
+ only at collocation or integration points.                           +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 50                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value
50
The solution is saved on a uniform grid of NX+1 points, covering the
interval (XA,XB).  Enter values for XA,XB.  These variables are usually
defaulted.

The defaults are XA = XGRID(1), XB = XGRID(NXGRID)

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default XA,XB                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
XA =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XB =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

The solution will be saved (for possible postprocessing) at the NSAVE+1
time points
T0 + K*(TF-T0)/NSAVE
K=0,...,NSAVE.  Enter a value for NSAVE.

If a user-specified constant time step is used, NSTEPS must be an
integer multiple of NSAVE.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NSAVE = 20                                                    \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NSAVE =
|---- Enter an integer value in the range 1 to +INFINITY
20
To generate tabular or graphical output, choose an output option from
the list below.

0. No further output is desired
1. Table of values at output points
The tabulated output is saved in a file.
2. One dimensional cross-sectional plots (versus X or T)
or, if applicable:
3. Surface plot of variable as function of X and T

Enter 0,1,2 or 3 to select an output option.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you decide later that you want additional types of plots not      +
+ requested during this interactive session, you will have to work     +
+ through a new interactive session, so it is recommended that you     +
+ request all output or plots you think you MIGHT eventually want now, +
+ during this session.                                                 +
+                                                                      +
+ Regardless of the options you select, a dummy subroutine POSTPR      +
+ will be included in the program created by the interactive driver;   +
+ you can add your own postprocessing code to this subroutine.         +
[RETURN]
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 3, the first time you see this message and                    \$
\$        0, the second time                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 0 to 3
3
Enter a value for IVAR, to select the variable to be plotted or
printed:
IVAR = 1 means H  (possibly as modified by UPRINT,..)
2       A
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ We want to plot H, so                                                \$
\$ enter: IVAR = 1                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
IVAR =
|---- Enter an integer value in the range 1 to 2
1
Enter the view latitude, VLAT, and the view longitude, VLON, desired
for this plot, in degrees.  VLAT and VLON must be between 10 and 80
degrees; each defaults to 45 degrees.  VLAT and VLON are usually
defaulted.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default VLAT and VLON                              \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
VLAT =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

VLON =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

Specify the range (UMIN,UMAX) for the dependent variable axis.  UMIN
and UMAX are often defaulted.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ By default, the plot will be scaled to just fit in the plot area.    +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default UMIN and UMAX                              \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
UMIN =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

UMAX =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

Enter a title, WITHOUT quotation marks.  A maximum of 40 characters
are allowed.  The default is no title.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: Water pressure head, H                                        \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
TITLE =             (Press [RETURN] to default)
|----Enter title or name---------------|
To generate tabular or graphical output, choose an output option from
the list below.

0. No further output is desired
1. Table of values at output points
The tabulated output is saved in a file.
2. One dimensional cross-sectional plots (versus X or T)
or, if applicable:
3. Surface plot of variable as function of X and T

Enter 0,1,2 or 3 to select an output option.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you decide later that you want additional types of plots not      +
+ requested during this interactive session, you will have to work     +
+ through a new interactive session, so it is recommended that you     +
+ request all output or plots you think you MIGHT eventually want now, +
+ during this session.                                                 +
+                                                                      +
+ Regardless of the options you select, a dummy subroutine POSTPR      +
+ will be included in the program created by the interactive driver;   +
+ you can add your own postprocessing code to this subroutine.         +
[RETURN]
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 3, the first time you see this message and                    \$
\$        0, the second time                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 0 to 3
0
Do you want to define any FORTRAN function subprograms used in any of
the FORTRAN 'expressions' entered earlier, entering them line by line?

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you selected double precision accuracy earlier, be sure to        +
+ declare these functions and their arguments DOUBLE PRECISION.        +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$ Then define the functions THETA and RK, one line at a time:          \$
\$       FUNCTION THETA(H,IDER)                                         \$
\$       IMPLICIT DOUBLE PRECISION (A-H,O-Z)                            \$
\$       B = 1.611 E6                                                   \$
\$       TS = 0.287                                                     \$
\$       TR = 0.075                                                     \$
\$       BETA = 3.96                                                    \$
\$       IF (H.LT.0.0) THEN                                             \$
\$          TH0 = TR + B*(TS-TR)/(B + ABS(H)**BETA)                     \$
\$          TH1 = B*(TS-TR)*BETA*ABS(H)**(BETA-1.)/(B+ABS(H)**BETA)**2  \$
\$       ELSE                                                           \$
\$          TH0 = TS                                                    \$
\$          TH1 = 0.0                                                   \$
[RETURN]
\$       ENDIF                                                          \$
\$       IF (IDER.EQ.0) THETA = TH0                                     \$
\$       IF (IDER.EQ.1) THETA = TH1                                     \$
\$       RETURN                                                         \$
\$       END                                                            \$
\$       FUNCTION RK(H)                                                 \$
\$       IMPLICIT DOUBLE PRECISION (A-H,O-Z)                            \$
\$       RKS = 0.00944                                                  \$
\$       A = 1.175 E6                                                   \$
\$       GAMMA = 4.74                                                   \$
\$       IF (H.LT.0.0) THEN                                             \$
\$          RK = RKS*A/(A+ABS(H)**GAMMA)                                \$
\$       ELSE                                                           \$
\$          RK = RKS                                                    \$
\$       ENDIF                                                          \$
\$       RETURN                                                         \$
\$       END                                                            \$
\$     [blank line]                                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
Remember to begin FORTRAN statements in column 7
|-----7-----Input FORTRAN now (type blank line to terminate)-----------|
FUNCTION THETA(H,IDER)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
B = 1.611 E6
TS = 0.287
TR = 0.075
BETA = 3.96
IF (H.LT.0.0) THEN
TH0 = TR + B*(TS-TR)/(B + ABS(H)**BETA)
TH1 = B*(TS-TR)*BETA*ABS(H)**(BETA-1.)/(B+ABS(H)**BETA)**2
ELSE
TH0 = TS
TH1 = 0.0
ENDIF
IF (IDER.EQ.0) THETA = TH0
IF (IDER.EQ.1) THETA = TH1
RETURN
END
FUNCTION RK(H)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
RKS = 0.00944
A = 1.175 E6
GAMMA = 4.74
IF (H.LT.0.0) THEN
RK = RKS*A/(A+ABS(H)**GAMMA)
ELSE
RK = RKS
ENDIF
RETURN
END

If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you want to define any FORTRAN function subprograms used in any   +
+ of the FORTRAN 'expressions' entered earlier, by interpolating the   +
+ tabular output saved in a file created on an earlier PDE2D run?      +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 14 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
More detailed information about PDE2D can be found in the book
"Solving Partial Differential Equation Applications with PDE2D,"
Granville Sewell, John Wiley & Sons, 2018.

*******************************************
***** Input program has been created  *****
*******************************************

```