Diffusion schemes -- grid centers


   subroutine diff_center(N,dt,cnpar,posconc,h,Bcup,Bcdw, &

This subroutine solves the one-dimensional diffusion equation including source terms,

$\displaystyle \dfrac{\partial {Y}}{\partial {t}} = \dfrac{\partial {}}{\partial...
... - \frac{1}{\tau_R}(Y-Y_{obs}) + Y L_{\text{sour}} + Q_{\text{sour}} \; , \quad$ (233)

for al variables defined at the centers of the grid cells, and a diffusion coefficient $ \nu_Y$ defined at the faces. Relaxation with time scale $ \tau_R$ towards observed values $ Y_{\text{obs}}$ is possible. $ L_{\text{sour}}$ specifies a linear source term, and $ Q_{\text{sour}}$ a constant source term. Central differences are used to discretize the problem as discussed in section 3.1.2. The diffusion term, the linear source term, and the linear part arising from the relaxation term are treated with an implicit method, whereas the constant source term is treated fully explicit.

The input parameters Bcup and Bcdw specify the type of the upper and lower boundary conditions, which can be either Dirichlet or Neumann-type. Bcup and Bcdw must have integer values corresponding to the parameters Dirichlet and Neumann defined in the module util, see section 8.2. Yup and Ydw are the values of the boundary conditions at the surface and the bottom. Depending on the values of Bcup and Bcdw, they represent either fluxes or prescribed values. The integer posconc indicates if a quantity is non-negative by definition (posconc=1, such as for concentrations) or not (posconc=0). For posconc=1 and negative boundary fluxes, the source term linearisation according to Patankar (1980) is applied.

Note that fluxes entering a boundary cell are counted positive by convention. The lower and upper position for prescribing these fluxes are located at the lowest und uppermost grid faces with index "0" and index "N", respectively. If values are prescribed, they are located at the centers with index "1" and index "N", respectively.


   use util,          only  : Dirichlet, Neumann
   use mtridiagonal
   number of vertical layers
   integer,  intent(in)                :: N
   time step (s)
   REALTYPE, intent(in)                :: dt
   "implicitness" parameter
   REALTYPE, intent(in)                :: cnpar
   1: non-negative concentration, 0: else
   integer, intent(in)                 :: posconc
   layer thickness (m)
   REALTYPE, intent(in)                :: h(0:N)
   type of upper BC
   integer,  intent(in)                :: Bcup
   type of lower BC
   integer,  intent(in)                :: Bcdw
   value of upper BC
   REALTYPE, intent(in)                :: Yup
   value of lower BC
   REALTYPE, intent(in)                :: Ydw
   diffusivity of Y
   REALTYPE, intent(in)                :: nuY(0:N)
   linear source term
   (treated implicitly)
   REALTYPE, intent(in)                :: Lsour(0:N)
   constant source term
   (treated explicitly)
   REALTYPE, intent(in)                :: Qsour(0:N)
   relaxation time (s)
   REALTYPE, intent(in)                :: Taur(0:N)
   observed value of Y
   REALTYPE, intent(in)                :: Yobs(0:N)
   REALTYPE                            :: Y(0:N)
   Original author(s): Lars Umlauf

Karsten Bolding 2012-12-28