c
c     file mud3cr.d
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 2008 by UCAR                 .
c  .                                                             .
c  .       UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                      MUDPACK version 5.0.1                  .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c
c ... file mud3cr.d
c
C     contains documentation for:
c     subroutine mud3cr(iparm,fparm,work,coef,bnd3cr,rhs,phi,mgopt,
c    +icros,crsxy,crsxz,crsyz,tol,maxit,iouter,rmax,ierror)
c     A sample fortran driver is file "tmud3cr.f".
c
c ... author and specialist
c
c          John C. Adams (National Center for Atmospheric Research)
c          email: mudpack.john@earthlink.net
c
c ... required MUDPACK files
c
c     mudcom.f
c
c
c ... purpose
c
c     subroutine mud3cr automatically discretizes and attempts to compute
c     the second order finite difference approximation to a three-
c     dimensional linear nonseparable elliptic partial differential
c     equation with cross derivative terms on a box.  the approximation
c     is generated on a uniform grid covering the box (see mesh description
c     below).  boundary conditions may be any combination of oblique mixed
c     derivative (see bnd3cr description below), specified (Dirchlet) or
c     periodic.  the form of the pde in operator notation is
c
c          l(p) + lxyz(p) = r(x,y,z)
c
c     where
c
c          l(p) = cxx(x,y,z)*pxx + cyy(x,y,z)*pyy + czz(z,y,z)*pzz +
c
c                 cx(x,y,z)*px + cy(x,y,z)*py + cz(x,y,z)*pz +
c
c                 ce(x,y,z)*p(x,y,z) = r(x,y,z)
c
c     and
c
c       lxyz(p) = cxy(x,y,z)*pxy + cxz(x,y,z)*pxz + cyz(x,y,z)*pyz
c
c     here cxx,cyy,czz,cx,cy,cz,ce,cxy,cxz,cyz  are the known real
c     coefficients of the pde; pxx,pyy,pzz,px,py,pz are the second and
c     first partial derivatives of the unknown solution function p(x,y,z)
c     with respect to the independent variables x,y,z; pxy,pxz, and pyz
c     are the second order mixed partial derivatives of p with respect
c     to xy,xz, and yz.  r(x,y,z) is the known right hand side of the pde.
c
c
c ... mesh description . . .
c
c     the approximation is generated on a uniform nx by ny by nz grid.
c     the grid is superimposed on the rectangular solution region
c
c          [xa,xb] x [yc,yd] x [ze,zf].
c
c     let
c
c          dlx = (xb-xa)/(nx-1), dly = (yd-yc)/(ny-1), dlz = (zf-ze)/(nz-1)
c
c     be the uniform grid increments in the x,y,z  directions. then
c
c          xi=xa+(i-1)*dlx,  yj=yc+(j-1)*dly, zk = ze+(k-1)*dlz
c
c     for i=1,...,nx; j=1,...,ny; k=1,...,nz  denote the x,y,z uniform
c     mesh points.
c
c
c
c ... methods
c
c
c     subroutine mud3cr is a recent addition to mudpack.  details
c     of the methods employeed by the other solvers in mudpack are in
c     [1,9,10].  [1,2,7,9,10] contain performance measurements on a variety
c     of elliptic pdes (see "references" in the file "readme").  the multi-
c     grid methods are described in documentation for the other solvers.
c *** mud3cr differs fundamentally from the other solvers in mudpack.
c     the full pde including cross derivative terms is discretized on
c     the INTERIOR of the solution region:
c
c       xa < x < xb, yc < y < yd, ze < z < zf
c
c     however, on nonspecified (nondirchlet) boundaries only l(p) is
c     discretized and the cross derivative term lxyz(p) is moved to the
c     right hand side of the pde and approximated by second order finite
c     finite difference formula applied to a previous estimate in p(k-1).
c     similarly, oblique mixed derivative boundary conditions (see bnd3cr)
c     are converted to a "mud3" type mixed normal form using second-order
c     finite difference formula applied to a previous estimate p(k-1) to
c     approximate non-normal derivative components.  for example if
c     the mixed derivative condition
c
c       py + a(x,z)*px + b(x,z)*pz + c(x,z)*p(x,yd,z) = gyd(y,z)
c
c     is specifed on the (x,z) plane of the upper y=yd boundary (see
c     bnd3cr for kbdy=4 below) then mud3cr converts this to the mixed
c     normal derivative form
c
c       py + c(x,z)*p(x,yd,z) = h(k,x,z)
c
c     where the modified right hand side h(k,x,z) is given by
c
c       h(k,x,z) = gyd(x,z) - [a(x,z)*dx(p(k-1)) + b(x,z)*dz(p(k-1)].
c
c     dx(p(k-1)) and dz(p(k-1)) are second order finite difference
c     approximations to the nonnormal partial derivatives px,pz using the
c     previous estimate in p(k-1).
c
c     the result of full discretization on interior grid points and partial
c     discretization with right hand side modifications on boundaries,
c     is a linear system which we denote by
c
c       D(p(k)) = r - Dxyz(p(k-1)).
c
c     D is the coefficient matrix coming from the discretization and
c     Dxyz(p(k-1)) stands for the right hand side modification obtained
c     by approximating boundary cross derivative terms and/or nonnormal
c     derivative components from mixed derivative boundary conditions
c     with second order finite difference formula applied to p(k-1).
c     with this notation, we formally describe the outer iteration employeed
c     by mud3cr:
c
c     algorithm mud3cr
c     .
c     set k = 0
c     .
c     set p(0) = 0.0 for all nonspecified grid points
c     .
c     repeat
c
c     .. k = k+1
c
c     .. solve D(p(k)) = r - Dxyz(p(k-1)) using multigrid iteration
c
c     .. set rmax(k) = ||p(k) - p(k-1)|| / ||p(k)||
c
c     until (rmax(k) < tol or k = maxit)
c     .
c     end mud3cr
c
c     tol is an error tolerance for convergence and maxit is a limit on
c     the number of outer iterations.  both are user prescribed input
c     arguments to mud3cr.  the maximum vector norm || || is used in
c     computing the relative difference between successive estimates in
c     rmax(k).  large values for maxit should not be used.
c
c *** note
c
c     originally a code, mud3cr0, was designed by moving all cross terms
c     to the right hand side and solving
c
c        l(p(k)) = r - lxyz(p(k-1)) for k=1,...
c
c     over the entire solution region including the interior.  in this
c     relatively straightforward approach, the standard mudpack solver
c     mud3 is used iteratively at each step.  however, convergence with
c     mud3cro is slow and unreliable.  an attempt was then made to
c     discretize the complete pde including cross terms and boundary
c     conditions over the entire region like the other solvers in mudpack.`
c     undoubtedly, this would have the most efficient and robust convergence
c     properties.  the main difficulty with this approach is the
c     unmanageable code complexity required to discretize all possible
c     combinations of nonzero cross derivative terms and oblique derivative
c     conditions.  for example, in the presence of nonzero cross terms,
c     cornors which are at the intersections of oblique derivative boundary
c     conditions become too complex (for this person) to discretize.
c     detection of combinations of oblique derivative conditions at cornors
c     which are singular (i.e., for which discretization leads to division
c     by 0.0) is a logical nightmare.
c
c     the present mud3cr is a middle ground between these two approaches.
c     it is an attempt to include as much of the pde cross terms as possible
c     in the discretization while bypassing the code complexity required to
c     include all possible boundary situations in the discretization.
c
c     by including the (manageable) discretization of cross terms on the
c     interior, the unfavorable convergence properties of the first approach
c     are (hopefully) avoided and the favorable convergence properties
c     of the second approach are (hopefully) obtained.  by moving nonzero
c     boundary cross terms and nonnormal components of mixed derivative
c     boundary conditions to the right hand side, the problems of too much
c     code complexity with discretization in the second approach are bypassed.
c
c     extensive numerical testing indicates for most problems mud3cr is more
c     robust than mud3cro.  convergence is more "iffy" and computationally
c     expensive than with the other solvers in mudpack.
c
c ... language                                                                  
c
c     fortran90/fortran77
c
c
c ... portability                                                               
c
c     mudpack5.0.1 software has been compiled and tested with Fortran77
c     and Fortran90 on a variety of platforms.
c                                                                               
c
c ... references (partial list)
c
c     for a complete list see "references" in the mudpack information and
c     directory file "readme"
c
c     [1] J. Adams, "MUDPACK: Multigrid Fortran Software for the Efficient
c     Solution of Linear Elliptic Partial Differential Equations,"
c     Applied Math. and Comput. vol.34, Nov 1989, pp.113-146.
c
c     [2] J. Adams,"FMG Results with the Multigrid Software Package MUDPACK,"
c     proceedings of the fourth Copper Mountain Conference on Multigrid, SIAM,
c     1989, pp.1-12.
c     .
c     .
c     .
c     [7] J. Adams, R. Garcia, B. Gross, J. Hack, D. Haidvogel, and V. Pizzo,
c     "Applications of Multigrid Software in the Atmospheric Sciences,"
c      Mon. Wea. Rev.,vol. 120 # 7, July 1992, pp. 1447-1458.
c     .
c     .
c     .
c     [9] J. Adams, "Recent Enhancements in MUDPACK, a Multigrid Software
c     package for Elliptic Partial Differential Equations," Applied Math.
c     and Comp., 1991, vol. 43, May 1991, pp. 79-94.
c
c     [10]J. Adams, "MUDPACK-2: Multigrid Software for Approximating
c     Elliptic Partial Differential Equations on Uniform Grids with
c     any Resolution," Applied Math. and Comp., 1993, vol. 53, February
c     1993, pp. 235-249
c     .
c     .
c     .
c
c **********************************************************************
c *** arguments ********************************************************
c **********************************************************************
c
c     arguments iparm,fparm,work,rhs,phi,coef,mgopt are the same as
c     those input to mud3 (see mud3.d for a detailed description) with the
c     following provisions:
c
c     (1) the minimum required work space length for mud3cr is increased
c         by approximately
c
c              nx*ny*nz*(1+8*(icros(1)+icros(2)+icros(3))/7 +
c
c              2*(icros(1)+icros(2)+icros(3))*(nx*ny+nx*nz+ny*nz)
c
c         words over the minimum work space required by mud3 (see icros
c         description below).  the exact minimal work space required
c         by mud3cr for the current set of input arguments is output
c         in iparm(22). * The exact minimal work length required
c         for the current method and grid size arguments can be
c         predetermined by calling mud3cr with iparm(21)=0 and
c         printout of iparm(22) or (in fortran 90 codes) dynamically
c         allocating work using the the value in iparm(22) in subsequent
c         calls to mud3cr.
c
c     (2) at least two calls to mud3cr are necessary to generate an
c         approximation. intl=iparm(1)=0 is required on the first
c         call.  this call will do "once only" discretization, and
c         set intermediate values in work which must be preserved
c         for noninitial calls.
c
c     (3) maxcy = iparm(18) must be 1 or 2 (see ierror = 13).
c
c     (4) tolmax = fparm(5) = 0.0 is required.  no "internal" error control
c         is allowed within multigrid cycling (see mud3.d)
c
c     (5) mgopt(1) = 0 is required.  only the default multigrid
c         options (W(2,1) cycles with cubic prolongation) can be used
c         with mud3cr
c
c *** new arguments
c
c      the arguments: bnd3cr,icros,crsxy,crsxz,crsyz,tol,maxit,iouter,rmax
c      are all new to mud3cr.  the error argument, ierror, has been expanded.
c      these are all described below:
c
c
c ... bnd3cr(kbdy,xory,yorz,a,b,c,g)
c
c      a subroutine with input arguments kbdy,xory,yorz and output
c      arguments a,b,c,g.  bnd3cr inputs OBLIQUE mixed derivative
c      conditions at any of the six x,y,z boundaries to mud3cr as
c      described below:
c
c         (1) the kbdy=1 boundary
c
c         this is the (y,z) plane x=xa where nxa=iparm(2)=2 flags
c         an oblique mixed boundary condition of the form
c
c           px + axa(y,z)*py + bxa(y,z)*pz +cxa(y,z)*p(xa,y,z) = gxa(y,z)
c
c         in this case kbdy=1,xory=y,yorz=z will be input to bnd3cr and
c         a,b,c,g corresponding to the known coefficients axa(y,z),bxa(y,z),
c         cxa(y,z),gxa(y,z) must be returned
c
c
c         (2) the kbdy=2 boundary
c
c         this is the (y,z) plane x=xb where nxb=iparm(3)=2 flags
c         an oblique mixed boundary condition of the form
c
c           px + axb(y,z)*py + bxb(y,z)*pz +cxb(y,z)*p(xb,y,z) = gxb(y,z)
c
c         in this case kbdy=2,xory=y,yorz=z will be input to bnd3cr and
c         a,b,c,g corresponding to the known coefficients axb(y,z),bxb(y,z),
c         cxb(y,z),gxb(y,z) must be returned
c
c         (3) the kbdy=3 boundary
c
c         this is the (x,z) plane y=yc where nyc=iparm(4)=2 flags
c         an oblique mixed boundary condition of the form
c
c           py + ayc(x,z)*px + byc(x,z)*pz +cyc(x,z)*p(x,yc,z) = gyc(x,z)
c
c         in this case kbdy=3,xory=x,yorz=z will be input to bnd3cr and
c         a,b,c,g corresponding to the known coefficients ayc(x,z),byc(x,z),
c         cyc(x,z),gyc(x,z) must be returned
c
c
c         (4) the kbdy=4 boundary
c
c         this is the (x,z) plane y=yd where nyd=iparm(5)=2 flags
c         an oblique mixed boundary condition of the form
c
c           py + ayd(x,z)*px  + byd(x,z)*pz +cyd(x,z)*p(x,yd,z) = gyd(x,z)
c
c         in this case kbdy=4,xory=x,yorz=z will be input to bnd3cr and
c         a,b,c,g corresponding to the known coefficients ayd(x,z),byd(x,z),
c         cyd(x,z),gyd(x,z) must be returned
c
c         (5) the kbdy=5 boundary
c
c         this is the (x,y) plane z=ze where nze=iparm(6)=2 flags
c         an oblique mixed boundary condition of the form
c
c           pz + aze(x,y)*px + bze(x,y)*py + cze(x,y)*p(x,y,ze) = gze(x,y)
c
c         in this case kbdy=5,xory=x,yorz=y will be input to bnd3cr and
c         a,b,c,g corresponding to the known coefficients aze(x,y),bze(x,y),
c         cze(x,y),gze(x,y) must be returned
c
c         (6) the kbdy=6 boundary
c
c         this is the (x,y) plane z=zf where nzf=iparm(7)=2 flags
c         an oblique mixed boundary condition of the form
c
c           pz + azf(x,y)*px + bzf(x,y)*py + czf(x,y)*p(x,y,zf) = gzf(x,y)
c
c         in this case kbdy=6,xory=x,yorz=y will be input to bnd3cr and
c         a,b,c,g corresponding to the known coefficients azf(x,y),bzf(x,y),
c         czf(x,y),gzf(x,y) must be returned
c
c
c     bnd3cr must be delcared "external" in the routine calling mud3cr
c     where its name may be different.  bnd3cr must be entered as a
c     dummy subroutine even if there are no derivative boundary conditions.
c     for an example of how to set up a subroutine to input derivative
c     boundary conditions, see the test program tmud3cr.f
c
c ... icros
c
c     an integer vector argument dimensioned 3 which flags the presence
c     or absence of cross derivative terms in the pde as follows:
c
c       icros(1) = 1 if cxy(x,y,z) is nonzero for any grid point (x,y,z)
c       icros(1) = 0 if cxy(x,y,z) = 0.0 for all grid points (x,y,z)
c
c       icros(2) = 1 if cxz(x,y,z) is nonzero for any grid point (x,y,z)
c       icros(2) = 0 if cxz(x,y,z) = 0.0 for all grid points (x,y,z)
c
c       icros(3) = 1 if cyz(x,y,z) is nonzero for any grid point (x,y,z)
c       icros(3) = 0 if cyz(x,y,z) = 0.0 for all grid points (x,y,z)
c
c
c ... crsxy(x,y,z,cxy)
c
c     if icros(1) = 1 then crsxy is a subroutine with arguments
c     (x,y,z,cxy) which supplies the xy cross derivative coefficient
c     cxy at the grid point (x,y,z).  if icros(1) = 0 then crsxy
c     is a dummy subroutine argument (i.e., it must be provided but
c     will not be invoked).
c
c
c ... crsxz(x,y,z,cxz)
c
c     if icros(2) = 1 then crsxz is a subroutine with arguments
c     (x,y,z,cxz) which supplies the xz cross derivative coefficient
c     cxz at the grid point (x,y,z).  if icros(2) = 0 then crsxz
c     is a dummy subroutine argument (i.e., it must be provided but
c     will not be invoked).
c
c
c ... crsyz(x,y,z,cyz)
c
c     if icros(3) = 1 then crsyz is a subroutine with arguments
c     (x,y,z,cyz) which supplies the yz cross derivative coefficient
c     cxy at the grid point (x,y,z).  if icros(3) = 0 then crsyz
c     is a dummy subroutine argument (i.e., it must be provided but
c     will not be invoked).
c
c     crsxy,crsxz,crsyz must be declared "external" in the routine
c     calling mud3cr.  the names chosen for these routines can be
c     different (see tmud3cr.f for an example)
c
c ... tol
c
c     tol is an error control argument for the outer iteration employed
c     by mud3cr (see "methods" description above).  if tol > 0.0 is input
c     then tol is a relative error tolerance for convergence.  the outer
c     iteration terminates and convergence is deemed to have occurred at the
c     k(th) iterate if the maximum relative difference, rmax(k), satisfies
c
c                  def
c           rmax(k) = ||p(k) - p(k-1)||/ ||p(k)|| < tol.
c
c     the last approximation p(maxit) is returned in phi even if
c     convergence does not occurr.  the maximum norm || || is used.
c     when tol = 0.0 is input, error control is not implemented and
c     exactly maxit (see below) outer iterations are executed in mud3cr.
c     the tol = 0.0 option eliminates unnecessary computation when
c     the user is certain of the required value for maxit.
c
c
c ... maxit
c
c     a limit on the outer iteration loop (see "method" description)
c     used to approximate the 3-d pde with cross derivative terms when
c     tol > 0.0.  if tol = 0.0 is entered then exactly maxit outer
c     iterations are performed and only rmax(maxit) is computed.  the
c     total number of relaxation sweeps performed at the finest grid
c     level is bounded by 3*maxcy*maxit.  large values for maxit should
c     not be used.
c
c
c ***********************************************************************
c ****output arguments**************************************************
c ***********************************************************************
c
c
c ... iparm(22)
c
c          on output iparm(22) contains the actual work space length
c          required by mud3cr for the current grid sizes and method.
c          this will be approximately

c              nx*ny*nz*(1+8*(icros(1)+icros(2)+icros(3))/7 +
c
c              2*(icros(1)+icros(2)+icros(3))*(nx*ny+nx*nz+ny*nz)
c
c          words longer than the space required by mud3 (see mud3.d)
c
c
c ... work
c
c          on output work contains intermediate values that must not be
c          destroyed if mud3cr is to be called again with iparm(1)=1
c          and iparm(17)=1.
c
c
c ... phi
c
c          on output phi(i,j,k) contains the approximation to
c          p(xi,yj,zk) for all mesh points i=1,...,nx; j=1,...,ny;
c          k=1,...,nz.  the last computed iterate in phi is returned
c          even if convergence is not obtained.
c
c
c ... iouter
c
c          the number of outer iterations (see "method" description above)
c          executed by mud3cr for the current call. maxit is an upper bound
c          for iouter
c
c
c ... rmax (see tol,maxit descriptions above)
c
c          a maxit dimensioned real vector.  if tol > 0.0 is input then
c          rmax(k) for k=1,...,iouter contain the maximum relative
c          difference between successive estimates.  rmax(k) is
c          given by
c
c               rmax(k) = ||p(k) - p(k-1)||/ ||p(k)||
c
c          for k=1,...,iouter.  the maximum norm || || is used.  either
c          iouter < maxit (convergence) or iouter = maxit is possible.
c          if tol = 0.0 input then exactly maxit outer iterations are
c          executed and only rmax(maxit) is computed.  in this case
c          rmax(1),...,rmax(maxit-1) are set to 0.0.  the tol = 0.0
c          option eliminates unnecessary computation when the user is
c          certain of the required value for maxit.
c
c
c ... ierror
c
c          an integer error argument which indicates fatal errors when
c          returned positive.  the negative values -5,-4,-3,-2,-1  and
c          ierror = 2,3,4,5,6,9,10 have the same meaning as described for
c          for mud3 (see mud3.d).  in addition:
c
c ***      new nonfatal error
c
c          ierror = -10 if tol > 0.0 is input (error control) and convergence
c          fails in maxit outer iterations.  in this case the latest
c          approximation p(maxit) is returned in phi (mud3cr can be recalled
c          with iparm(1)=iparm(17)=1 to improve the approximation as long
c          as all other arguments are unchanged)
c
c ***      new fatal errors
c
c ... ierror
c
c     = 1 if intl=iparm(1) is not 0 on initial call or not 0 or 1
c         on subsequent calls of if intl=0 and iguess=iparm(17)=1
c
c     = 7 if maxcy = iparm(18) is not 1 or 2
c
c     = 8 if method = iparm(19) is less than 0 or greater than 7
c         mud3cr does not allow planar relaxation.  meth2=iparm(20)
c         is not used or checked.
c
c     =11 if tolmax = fparm(7) is not 0.0
c
c     =12 if kcycle = mgopt(1) is not 0
c
c     =13 if icros(1) or icros(2) or icros(3) is not 0 or 1
c
c     =14 if tol < 0.0
c
c     =15 if maxit < 1
c
c  ***********************************************************************
c  ***********************************************************************
c
c     end of mud3cr documentation
c
c  ***********************************************************************
c  ***********************************************************************
c