/*****************************
 *  FUNCTION:  load_AF_CF_J  *
 *****************************/
void load_AF_CF_J (F, j, unkSW, unkSE, ccx, ccy, cx, cy, dx, dy)
struct triblock22 *F;
int j, unkSW, unkSE;
double **ccx, **ccy, **cx, **cy, *dx, *dy;
{
int i;

if (unkSW == 0)
    {
    F->a[0].el[0][0] = J_fRx_fRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[0]);
    F->a[0].el[1][0] = J_fRx_fRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[0]);
    }
else if (unkSW == 1)
    {
    F->a[0].el[0][0] = J_gRx_fRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[0]);
    F->a[0].el[1][0] = J_gRx_fRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[0]);
    }
else if (unkSW == 2)
    {
    F->a[0].el[0][0] = J_fRx_gRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[0]);
    F->a[0].el[1][0] = J_fRx_gRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[0]);
    }
else
    {
    F->a[0].el[0][0] = J_gRx_gRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[0]);
    F->a[0].el[1][0] = J_gRx_gRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[0]);
    }

if (unkSE == 0)
    {
    F->a[xm-1].el[0][1] = J_fLx_fRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[0]);
    F->a[xm-1].el[1][1] = J_fLx_fRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[0]);
    }
else if (unkSE == 1)
    {
    F->a[xm-1].el[0][1] = J_gLx_fRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[0]);
    F->a[xm-1].el[1][1] = J_gLx_fRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[0]);
    }
else if (unkSE == 2)
    {
    F->a[xm-1].el[0][1] = J_fLx_gRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[0]);
    F->a[xm-1].el[1][1] = J_fLx_gRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[0]);
    }
else
    {
    F->a[xm-1].el[0][1] = J_gLx_gRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[0]);
    F->a[xm-1].el[1][1] = J_gLx_gRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[0]);
    }
  
for (i = 0; i <= xm-2; i++)
    {
    if (BCTypeS == 1)
        {
	F->a[i]  .el[0][1] = J_fLx_gRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[0]);
	F->a[i]  .el[1][1] = J_fLx_gRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[0]);
	F->b[i]  .el[0]    = J_gLx_gRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[0]);
	F->b[i]  .el[1]    = J_gLx_gRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[0]);
	F->c[i]  .el[0]    = J_fRx_gRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[0]);
	F->c[i]  .el[1]    = J_fRx_gRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[0]);
	F->a[i+1].el[0][0] = J_gRx_gRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[0]);
	F->a[i+1].el[1][0] = J_gRx_gRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[0]);
        }
    else
        {
	F->a[i]  .el[0][1] = J_fLx_fRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[0]);
	F->a[i]  .el[1][1] = J_fLx_fRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[0]);
	F->b[i]  .el[0]    = J_gLx_fRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[0]);
	F->b[i]  .el[1]    = J_gLx_fRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[0]);
	F->c[i]  .el[0]    = J_fRx_fRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[0]);
	F->c[i]  .el[1]    = J_fRx_fRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[0]);
	F->a[i+1].el[0][0] = J_gRx_fRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[0]);
	F->a[i+1].el[1][0] = J_gRx_fRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[0]);
        }
    }
}  /* end of function load_AF_CF_J */


/*****************************
 *  FUNCTION:  load_AF_CF_L  *
 *****************************/
void load_AF_CF_L(F, j, unkSW, unkSE, ccx, ccy, cx, cy, dx, dy, t)
struct triblock22 *F;
int j, unkSW, unkSE;
double **ccx, **ccy, **cx, **cy, *dx, *dy, t;
{
int i;

if (unkSW == 0)
    {
    F->a[0].el[0][0] = L_fRx_fRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[0] ,t);
    F->a[0].el[1][0] = L_fRx_fRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[0] ,t);
    }
else if (unkSW == 1)
    {
    F->a[0].el[0][0] = L_gRx_fRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[0] ,t);
    F->a[0].el[1][0] = L_gRx_fRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[0] ,t);
    }
else if (unkSW == 2)
    {
    F->a[0].el[0][0] = L_fRx_gRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[0] ,t);
    F->a[0].el[1][0] = L_fRx_gRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[0] ,t);
    }
else
    {
    F->a[0].el[0][0] = L_gRx_gRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[0] ,t);
    F->a[0].el[1][0] = L_gRx_gRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[0] ,t);
    }

if (unkSE == 0)
    {
    F->a[xm-1].el[0][1] = L_fLx_fRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[0] ,t);
    F->a[xm-1].el[1][1] = L_fLx_fRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[0] ,t);
    }
else if (unkSE == 1)
    {
    F->a[xm-1].el[0][1] = L_gLx_fRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[0] ,t);
    F->a[xm-1].el[1][1] = L_gLx_fRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[0] ,t);
    }
else if (unkSE == 2)
    {
    F->a[xm-1].el[0][1] = L_fLx_gRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[0] ,t);
    F->a[xm-1].el[1][1] = L_fLx_gRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[0] ,t);
    }
else
    {
    F->a[xm-1].el[0][1] = L_gLx_gRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[0] ,t);
    F->a[xm-1].el[1][1] = L_gLx_gRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[0] ,t);
    }
  
for (i = 0; i <= xm-2; i++)
    {
    if (BCTypeS == 1)
        {
	F->a[i]  .el[0][1] = L_fLx_gRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[0] ,t);
	F->a[i]  .el[1][1] = L_fLx_gRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[0] ,t);
	F->b[i]  .el[0]    = L_gLx_gRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[0] ,t);
	F->b[i]  .el[1]    = L_gLx_gRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[0] ,t);
	F->c[i]  .el[0]    = L_fRx_gRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[0] ,t);
	F->c[i]  .el[1]    = L_fRx_gRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[0] ,t);
	F->a[i+1].el[0][0] = L_gRx_gRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[0] ,t);
	F->a[i+1].el[1][0] = L_gRx_gRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[0] ,t);
        }
    else
        {
	F->a[i]  .el[0][1] = L_fLx_fRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[0] ,t);
	F->a[i]  .el[1][1] = L_fLx_fRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[0] ,t);
	F->b[i]  .el[0]    = L_gLx_fRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[0] ,t);
	F->b[i]  .el[1]    = L_gLx_fRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[0] ,t);
	F->c[i]  .el[0]    = L_fRx_fRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[0] ,t);
	F->c[i]  .el[1]    = L_fRx_fRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[0] ,t);
	F->a[i+1].el[0][0] = L_gRx_fRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[0] ,t);
	F->a[i+1].el[1][0] = L_gRx_fRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[0] ,t);
        }
    }
}  /* end of function load_AF_CF_L */


/*****************************
 *  FUNCTION:  load_AL_BL_J  *
 *****************************/
void load_AL_BL_J(L, j, unkNW, unkNE, ccx, ccy, cx, cy, dx, dy)
struct triblock22 *L;
int j, unkNW, unkNE;
double **ccx, **ccy, **cx, **cy, *dx, *dy;
{
int i;

if (unkNW == 0)
    {
    L->a[0].el[0][0] = J_fRx_fLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[ym-1]);
    L->a[0].el[1][0] = J_fRx_fLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[ym-1]);
    }
else if (unkNW == 1)
    {
    L->a[0].el[0][0] = J_gRx_fLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[ym-1]);
    L->a[0].el[1][0] = J_gRx_fLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[ym-1]);
    }
else if (unkNW == 2)
    {
    L->a[0].el[0][0] = J_fRx_gLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[ym-1]);
    L->a[0].el[1][0] = J_fRx_gLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[ym-1]);
    }
else
    {
    L->a[0].el[0][0] = J_gRx_gLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[ym-1]);
    L->a[0].el[1][0] = J_gRx_gLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[ym-1]);
    }

if (unkNE == 0)
    {
    L->a[xm-1].el[0][1] = J_fLx_fLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[ym-1]);
    L->a[xm-1].el[1][1] = J_fLx_fLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[ym-1]);
    }
else if (unkNE == 1)
    {
    L->a[xm-1].el[0][1] = J_gLx_fLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[ym-1]);
    L->a[xm-1].el[1][1] = J_gLx_fLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[ym-1]);
    }
else if (unkNE == 2)
    {
    L->a[xm-1].el[0][1] = J_fLx_gLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[ym-1]);
    L->a[xm-1].el[1][1] = J_fLx_gLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[ym-1]);
    }
else
    {
    L->a[xm-1].el[0][1] = J_gLx_gLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[ym-1]);
    L->a[xm-1].el[1][1] = J_gLx_gLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[ym-1]);
    }

for (i = 0; i <= xm-2; i++)
    {
    if (BCTypeN == 1)
        {
	L->a[i]  .el[0][1] = J_fLx_gLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[ym-1]);
	L->a[i]  .el[1][1] = J_fLx_gLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[ym-1]);
	L->b[i]  .el[0]    = J_gLx_gLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[ym-1]);
	L->b[i]  .el[1]    = J_gLx_gLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[ym-1]);
	L->c[i]  .el[0]    = J_fRx_gLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[ym-1]);
	L->c[i]  .el[1]    = J_fRx_gLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[ym-1]);
	L->a[i+1].el[0][0] = J_gRx_gLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[ym-1]);
	L->a[i+1].el[1][0] = J_gRx_gLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[ym-1]);
        }
    else
        {
	L->a[i]  .el[0][1] = J_fLx_fLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[ym-1]);
	L->a[i]  .el[1][1] = J_fLx_fLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[ym-1]);
	L->b[i]  .el[0]    = J_gLx_fLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[ym-1]);
	L->b[i]  .el[1]    = J_gLx_fLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[ym-1]);
	L->c[i]  .el[0]    = J_fRx_fLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[ym-1]);
	L->c[i]  .el[1]    = J_fRx_fLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[ym-1]);
	L->a[i+1].el[0][0] = J_gRx_fLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[ym-1]);
	L->a[i+1].el[1][0] = J_gRx_fLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[ym-1]);
        }
    }
}  /* end of function load_AL_BL_J */


/*****************************
 *  FUNCTION:  load_AL_BL_L  *
 *****************************/
void load_AL_BL_L(L, j, unkNW, unkNE, ccx, ccy, cx, cy, dx, dy, t)
struct triblock22 *L;
int j, unkNW, unkNE;
double **ccx, **ccy, **cx, **cy, *dx, *dy, t;
{
int i;

if (unkNW == 0)
    {
    L->a[0].el[0][0] = L_fRx_fLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[ym-1] ,t);
    L->a[0].el[1][0] = L_fRx_fLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[ym-1] ,t);
    }
else if (unkNW == 1)
    {
    L->a[0].el[0][0] = L_gRx_fLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[ym-1] ,t);
    L->a[0].el[1][0] = L_gRx_fLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[ym-1] ,t);
    }
else if (unkNW == 2)
    {
    L->a[0].el[0][0] = L_fRx_gLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[ym-1] ,t);
    L->a[0].el[1][0] = L_fRx_gLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[ym-1] ,t);
    }
else
    {
    L->a[0].el[0][0] = L_gRx_gLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[ym-1] ,t);
    L->a[0].el[1][0] = L_gRx_gLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[ym-1] ,t);
    }

if (unkNE == 0)
    {
    L->a[xm-1].el[0][1] = L_fLx_fLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[ym-1] ,t);
    L->a[xm-1].el[1][1] = L_fLx_fLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[ym-1] ,t);
    }
else if (unkNE == 1)
    {
    L->a[xm-1].el[0][1] = L_gLx_fLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[ym-1] ,t);
    L->a[xm-1].el[1][1] = L_gLx_fLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[ym-1] ,t);
    }
else if (unkNE == 2)
    {
    L->a[xm-1].el[0][1] = L_fLx_gLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[ym-1] ,t);
    L->a[xm-1].el[1][1] = L_fLx_gLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[ym-1] ,t);
    }
else
    {
    L->a[xm-1].el[0][1] = L_gLx_gLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[ym-1] ,t);
    L->a[xm-1].el[1][1] = L_gLx_gLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[ym-1] ,t);
    }

for (i = 0; i <= xm-2; i++)
    {
    if (BCTypeN == 1)
        {
	L->a[i]  .el[0][1] = L_fLx_gLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[ym-1] ,t);
	L->a[i]  .el[1][1] = L_fLx_gLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[ym-1] ,t);
	L->b[i]  .el[0]    = L_gLx_gLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[ym-1] ,t);
	L->b[i]  .el[1]    = L_gLx_gLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[ym-1] ,t);
	L->c[i]  .el[0]    = L_fRx_gLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[ym-1] ,t);
	L->c[i]  .el[1]    = L_fRx_gLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[ym-1] ,t);
	L->a[i+1].el[0][0] = L_gRx_gLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[ym-1] ,t);
	L->a[i+1].el[1][0] = L_gRx_gLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[ym-1] ,t);
        }
    else
        {
	L->a[i]  .el[0][1] = L_fLx_fLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[ym-1] ,t);
	L->a[i]  .el[1][1] = L_fLx_fLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[ym-1] ,t);
	L->b[i]  .el[0]    = L_gLx_fLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[ym-1] ,t);
	L->b[i]  .el[1]    = L_gLx_fLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[ym-1] ,t);
	L->c[i]  .el[0]    = L_fRx_fLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[ym-1] ,t);
	L->c[i]  .el[1]    = L_fRx_fLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[ym-1] ,t);
	L->a[i+1].el[0][0] = L_gRx_fLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[ym-1] ,t);
	L->a[i+1].el[1][0] = L_gRx_fLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[ym-1] ,t);
        }
    }
}  /* end of function load_AL_BL_L */


/*************************
 *  FUNCTION:  load_B_J  *
 *************************/
load_B_J(B, j, k, ccx, ccy, cx, cy, dx, dy)
struct triblock24 *B;
int j, k;
double **ccx, **ccy, **cx, **cy, *dx, *dy;
{
int i;

if (BCTypeW == 1)
    {
    B->a[0].el[0].el[0][0] = J_gRx_fLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k]);
    B->a[0].el[0].el[0][1] = J_gRx_gLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k]);
    B->a[0].el[0].el[1][0] = J_gRx_fLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k]);
    B->a[0].el[0].el[1][1] = J_gRx_gLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k]);
    }
else
    {
    B->a[0].el[0].el[0][0] = J_fRx_fLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k]);
    B->a[0].el[0].el[0][1] = J_fRx_gLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k]);
    B->a[0].el[0].el[1][0] = J_fRx_fLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k]);
    B->a[0].el[0].el[1][1] = J_fRx_gLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k]);
    }

if (BCTypeE == 1)
    {
    B->a[xm-1].el[1].el[0][0] = J_gLx_fLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k]);
    B->a[xm-1].el[1].el[0][1] = J_gLx_gLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k]);
    B->a[xm-1].el[1].el[1][0] = J_gLx_fLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k]);
    B->a[xm-1].el[1].el[1][1] = J_gLx_gLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k]);
    }
else
    {
    B->a[xm-1].el[1].el[0][0] = J_fLx_fLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k]);
    B->a[xm-1].el[1].el[0][1] = J_fLx_gLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k]);
    B->a[xm-1].el[1].el[1][0] = J_fLx_fLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k]);
    B->a[xm-1].el[1].el[1][1] = J_fLx_gLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k]);
    }

for (i = 0; i <= xm-2; i++)
    {
    B->a[i]  .el[1].el[0][0] = J_fLx_fLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k]);
    B->a[i]  .el[1].el[0][1] = J_fLx_gLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k]);
    B->a[i]  .el[1].el[1][0] = J_fLx_fLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k]);
    B->a[i]  .el[1].el[1][1] = J_fLx_gLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k]);
    
    B->b[i]        .el[0][0] = J_gLx_fLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k]);
    B->b[i]        .el[0][1] = J_gLx_gLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k]);
    B->b[i]        .el[1][0] = J_gLx_fLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k]);
    B->b[i]        .el[1][1] = J_gLx_gLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k]);
    
    B->c[i]        .el[0][0] = J_fRx_fLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k]);
    B->c[i]        .el[0][1] = J_fRx_gLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k]);
    B->c[i]        .el[1][0] = J_fRx_fLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k]);
    B->c[i]        .el[1][1] = J_fRx_gLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k]);
    
    B->a[i+1].el[0].el[0][0] = J_gRx_fLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k]);
    B->a[i+1].el[0].el[0][1] = J_gRx_gLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k]);
    B->a[i+1].el[0].el[1][0] = J_gRx_fLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k]);
    B->a[i+1].el[0].el[1][1] = J_gRx_gLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k]);
    }
}  /* end of function load_B_J */


/*************************
 *  FUNCTION:  load_B_L  *
 *************************/
load_B_L(B, j, k, ccx, ccy, cx, cy, dx, dy, t)
struct triblock24 *B;
int j, k;
double **ccx, **ccy, **cx, **cy, *dx, *dy, t;
{
int i;

if (BCTypeW == 1)
    {
    B->a[0].el[0].el[0][0] = L_gRx_fLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k] ,t);
    B->a[0].el[0].el[0][1] = L_gRx_gLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k] ,t);
    B->a[0].el[0].el[1][0] = L_gRx_fLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k] ,t);
    B->a[0].el[0].el[1][1] = L_gRx_gLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k] ,t);
    }
else
    {
    B->a[0].el[0].el[0][0] = L_fRx_fLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k] ,t);
    B->a[0].el[0].el[0][1] = L_fRx_gLy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k] ,t);
    B->a[0].el[0].el[1][0] = L_fRx_fLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k] ,t);
    B->a[0].el[0].el[1][1] = L_fRx_gLy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k] ,t);
    }

if (BCTypeE == 1)
    {
    B->a[xm-1].el[1].el[0][0] = L_gLx_fLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k] ,t);
    B->a[xm-1].el[1].el[0][1] = L_gLx_gLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k] ,t);
    B->a[xm-1].el[1].el[1][0] = L_gLx_fLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k] ,t);
    B->a[xm-1].el[1].el[1][1] = L_gLx_gLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k] ,t);
    }
else
    {
    B->a[xm-1].el[1].el[0][0] = L_fLx_fLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k] ,t);
    B->a[xm-1].el[1].el[0][1] = L_fLx_gLy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k] ,t);
    B->a[xm-1].el[1].el[1][0] = L_fLx_fLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k] ,t);
    B->a[xm-1].el[1].el[1][1] = L_fLx_gLy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k] ,t);
    }

for (i = 0; i <= xm-2; i++)
    {
    B->a[i]  .el[1].el[0][0] = L_fLx_fLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k] ,t);
    B->a[i]  .el[1].el[0][1] = L_fLx_gLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k] ,t);
    B->a[i]  .el[1].el[1][0] = L_fLx_fLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k] ,t);
    B->a[i]  .el[1].el[1][1] = L_fLx_gLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k] ,t);
    
    B->b[i]        .el[0][0] = L_gLx_fLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k] ,t);
    B->b[i]        .el[0][1] = L_gLx_gLy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k] ,t);
    B->b[i]        .el[1][0] = L_gLx_fLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k] ,t);
    B->b[i]        .el[1][1] = L_gLx_gLy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k] ,t);
    
    B->c[i]        .el[0][0] = L_fRx_fLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k] ,t);
    B->c[i]        .el[0][1] = L_fRx_gLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k] ,t);
    B->c[i]        .el[1][0] = L_fRx_fLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k] ,t);
    B->c[i]        .el[1][1] = L_fRx_gLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k] ,t);
    
    B->a[i+1].el[0].el[0][0] = L_gRx_fLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k] ,t);
    B->a[i+1].el[0].el[0][1] = L_gRx_gLy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k] ,t);
    B->a[i+1].el[0].el[1][0] = L_gRx_fLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k] ,t);
    B->a[i+1].el[0].el[1][1] = L_gRx_gLy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k] ,t);
    }
}  /* end of function load_B_L */


/*************************
 *  FUNCTION:  load_C_J  *
 *************************/
load_C_J(C, j, k, ccx, ccy, cx, cy, dx, dy)
struct triblock24 *C;
int j, k;
double **ccx, **ccy, **cx, **cy, *dx, *dy;
{
int i;

if (BCTypeW == 1)
    {
    C->a[0].el[0].el[0][0] = J_gRx_fRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k]);
    C->a[0].el[0].el[0][1] = J_gRx_gRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k]);
    C->a[0].el[0].el[1][0] = J_gRx_fRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k]);
    C->a[0].el[0].el[1][1] = J_gRx_gRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k]);
    }
else
    {
    C->a[0].el[0].el[0][0] = J_fRx_fRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k]);
    C->a[0].el[0].el[0][1] = J_fRx_gRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k]);
    C->a[0].el[0].el[1][0] = J_fRx_fRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k]);
    C->a[0].el[0].el[1][1] = J_fRx_gRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k]);
    }

if (BCTypeE == 1)
    {
    C->a[xm-1].el[1].el[0][0] = J_gLx_fRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k]);
    C->a[xm-1].el[1].el[0][1] = J_gLx_gRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k]);
    C->a[xm-1].el[1].el[1][0] = J_gLx_fRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k]);
    C->a[xm-1].el[1].el[1][1] = J_gLx_gRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k]);
    }
else
    {
    C->a[xm-1].el[1].el[0][0] = J_fLx_fRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k]);
    C->a[xm-1].el[1].el[0][1] = J_fLx_gRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k]);
    C->a[xm-1].el[1].el[1][0] = J_fLx_fRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k]);
    C->a[xm-1].el[1].el[1][1] = J_fLx_gRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k]);
    }

for (i = 0; i <= xm-2; i++)
    {
    C->a[i]  .el[1].el[0][0] = J_fLx_fRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k]);
    C->a[i]  .el[1].el[0][1] = J_fLx_gRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k]);
    C->a[i]  .el[1].el[1][0] = J_fLx_fRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k]);
    C->a[i]  .el[1].el[1][1] = J_fLx_gRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k]);
    
    C->b[i]        .el[0][0] = J_gLx_fRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k]);
    C->b[i]        .el[0][1] = J_gLx_gRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k]);
    C->b[i]        .el[1][0] = J_gLx_fRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k]);
    C->b[i]        .el[1][1] = J_gLx_gRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k]);
    
    C->c[i]        .el[0][0] = J_fRx_fRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k]);
    C->c[i]        .el[0][1] = J_fRx_gRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k]);
    C->c[i]        .el[1][0] = J_fRx_fRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k]);
    C->c[i]        .el[1][1] = J_fRx_gRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k]);
    
    C->a[i+1].el[0].el[0][0] = J_gRx_fRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k]);
    C->a[i+1].el[0].el[0][1] = J_gRx_gRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k]);
    C->a[i+1].el[0].el[1][0] = J_gRx_fRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k]);
    C->a[i+1].el[0].el[1][1] = J_gRx_gRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k]);
    }
}  /* end of function load_C_J */


/*************************
 *  FUNCTION:  load_C_L  *
 *************************/
load_C_L(C, j, k, ccx, ccy, cx, cy, dx, dy, t)
struct triblock24 *C;
int j, k;
double **ccx, **ccy, **cx, **cy, *dx, *dy, t;
{
int i;

if (BCTypeW == 1)
    {
    C->a[0].el[0].el[0][0] = L_gRx_fRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k] ,t);
    C->a[0].el[0].el[0][1] = L_gRx_gRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k] ,t);
    C->a[0].el[0].el[1][0] = L_gRx_fRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k] ,t);
    C->a[0].el[0].el[1][1] = L_gRx_gRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k] ,t);
    }
else
    {
    C->a[0].el[0].el[0][0] = L_fRx_fRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k] ,t);
    C->a[0].el[0].el[0][1] = L_fRx_gRy(ccx[0][j], ccy[0][j], cx[0][j], cy[0][j], dx[0], dy[k] ,t);
    C->a[0].el[0].el[1][0] = L_fRx_fRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k] ,t);
    C->a[0].el[0].el[1][1] = L_fRx_gRy(ccx[1][j], ccy[1][j], cx[1][j], cy[1][j], dx[0], dy[k] ,t);
    }

if (BCTypeE == 1)
    {
    C->a[xm-1].el[1].el[0][0] = L_gLx_fRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k] ,t);
    C->a[xm-1].el[1].el[0][1] = L_gLx_gRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k] ,t);
    C->a[xm-1].el[1].el[1][0] = L_gLx_fRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k] ,t);
    C->a[xm-1].el[1].el[1][1] = L_gLx_gRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k] ,t);
    }
else
    {
    C->a[xm-1].el[1].el[0][0] = L_fLx_fRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k] ,t);
    C->a[xm-1].el[1].el[0][1] = L_fLx_gRy(ccx[xc-2][j], ccy[xc-2][j], cx[xc-2][j], cy[xc-2][j], dx[xm-1], dy[k] ,t);
    C->a[xm-1].el[1].el[1][0] = L_fLx_fRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k] ,t);
    C->a[xm-1].el[1].el[1][1] = L_fLx_gRy(ccx[xc-1][j], ccy[xc-1][j], cx[xc-1][j], cy[xc-1][j], dx[xm-1], dy[k] ,t);
    }

for (i = 0; i <= xm-2; i++)
    {
    C->a[i]  .el[1].el[0][0] = L_fLx_fRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k] ,t);
    C->a[i]  .el[1].el[0][1] = L_fLx_gRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k] ,t);
    C->a[i]  .el[1].el[1][0] = L_fLx_fRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k] ,t);
    C->a[i]  .el[1].el[1][1] = L_fLx_gRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k] ,t);
    
    C->b[i]        .el[0][0] = L_gLx_fRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k] ,t);
    C->b[i]        .el[0][1] = L_gLx_gRy(ccx[2*i]  [j], ccy[2*i]  [j], cx[2*i]  [j], cy[2*i]  [j], dx[i],   dy[k] ,t);
    C->b[i]        .el[1][0] = L_gLx_fRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k] ,t);
    C->b[i]        .el[1][1] = L_gLx_gRy(ccx[2*i+1][j], ccy[2*i+1][j], cx[2*i+1][j], cy[2*i+1][j], dx[i],   dy[k] ,t);
    
    C->c[i]        .el[0][0] = L_fRx_fRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k] ,t);
    C->c[i]        .el[0][1] = L_fRx_gRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k] ,t);
    C->c[i]        .el[1][0] = L_fRx_fRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k] ,t);
    C->c[i]        .el[1][1] = L_fRx_gRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k] ,t);
    
    C->a[i+1].el[0].el[0][0] = L_gRx_fRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k] ,t);
    C->a[i+1].el[0].el[0][1] = L_gRx_gRy(ccx[2*i+2][j], ccy[2*i+2][j], cx[2*i+2][j], cy[2*i+2][j], dx[i+1], dy[k] ,t);
    C->a[i+1].el[0].el[1][0] = L_gRx_fRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k] ,t);
    C->a[i+1].el[0].el[1][1] = L_gRx_gRy(ccx[2*i+3][j], ccy[2*i+3][j], cx[2*i+3][j], cy[2*i+3][j], dx[i+1], dy[k] ,t);
    }
}  /* end of function load_C_L */


/*************************
 *  FUNCTION:  load_A_J  *
 *************************/
load_A_J(A, j, k, ccx, ccy, cx, cy, dx, dy)
struct triblock44 *A;
int j, k;
double **ccx, **ccy, **cx, **cy, *dx, *dy;
{
int i;

if (BCTypeW == 1)
    {
    A->a[0].A.el[0][0] = J_gRx_fLy(ccx[0][j-1], ccy[0][j-1], cx[0][j-1], cy[0][j-1], dx[0], dy[k-1]);
    A->a[0].A.el[0][1] = J_gRx_gLy(ccx[0][j-1], ccy[0][j-1], cx[0][j-1], cy[0][j-1], dx[0], dy[k-1]);
    A->a[0].A.el[1][0] = J_gRx_fRy(ccx[0][j],   ccy[0][j],   cx[0][j],   cy[0][j]  , dx[0], dy[k]);
    A->a[0].A.el[1][1] = J_gRx_gRy(ccx[0][j],   ccy[0][j],   cx[0][j],   cy[0][j]  , dx[0], dy[k]);
    
    A->a[0].C.el[0][0] = J_gRx_fLy(ccx[1][j-1], ccy[1][j-1], cx[1][j-1], cy[1][j-1], dx[0], dy[k-1]);
    A->a[0].C.el[0][1] = J_gRx_gLy(ccx[1][j-1], ccy[1][j-1], cx[1][j-1], cy[1][j-1], dx[0], dy[k-1]);
    A->a[0].C.el[1][0] = J_gRx_fRy(ccx[1][j],   ccy[1][j],   cx[1][j],   cy[1][j]  , dx[0], dy[k]);
    A->a[0].C.el[1][1] = J_gRx_gRy(ccx[1][j],   ccy[1][j],   cx[1][j],   cy[1][j]  , dx[0], dy[k]);
    }
else
    {
    A->a[0].A.el[0][0] = J_fRx_fLy(ccx[0][j-1], ccy[0][j-1], cx[0][j-1], cy[0][j-1], dx[0], dy[k-1]);
    A->a[0].A.el[0][1] = J_fRx_gLy(ccx[0][j-1], ccy[0][j-1], cx[0][j-1], cy[0][j-1], dx[0], dy[k-1]);
    A->a[0].A.el[1][0] = J_fRx_fRy(ccx[0][j],   ccy[0][j],   cx[0][j],   cy[0][j]  , dx[0], dy[k]);
    A->a[0].A.el[1][1] = J_fRx_gRy(ccx[0][j],   ccy[0][j],   cx[0][j],   cy[0][j]  , dx[0], dy[k]);
    
    A->a[0].C.el[0][0] = J_fRx_fLy(ccx[1][j-1], ccy[1][j-1], cx[1][j-1], cy[1][j-1], dx[0], dy[k-1]);
    A->a[0].C.el[0][1] = J_fRx_gLy(ccx[1][j-1], ccy[1][j-1], cx[1][j-1], cy[1][j-1], dx[0], dy[k-1]);
    A->a[0].C.el[1][0] = J_fRx_fRy(ccx[1][j],   ccy[1][j],   cx[1][j],   cy[1][j]  , dx[0], dy[k]);
    A->a[0].C.el[1][1] = J_fRx_gRy(ccx[1][j],   ccy[1][j],   cx[1][j],   cy[1][j]  , dx[0], dy[k]);
    }

if (BCTypeE == 1)
    {
    A->a[xm-1].B.el[0][0] = J_gLx_fLy(ccx[xc-2][j-1], ccy[xc-2][j-1], cx[xc-2][j-1], cy[xc-2][j-1], dx[xm-1], dy[k-1]);
    A->a[xm-1].B.el[0][1] = J_gLx_gLy(ccx[xc-2][j-1], ccy[xc-2][j-1], cx[xc-2][j-1], cy[xc-2][j-1], dx[xm-1], dy[k-1]);
    A->a[xm-1].B.el[1][0] = J_gLx_fRy(ccx[xc-2][j],   ccy[xc-2][j],   cx[xc-2][j],   cy[xc-2][j],   dx[xm-1], dy[k]);
    A->a[xm-1].B.el[1][1] = J_gLx_gRy(ccx[xc-2][j],   ccy[xc-2][j],   cx[xc-2][j],   cy[xc-2][j],   dx[xm-1], dy[k]);
    
    A->a[xm-1].D.el[0][0] = J_gLx_fLy(ccx[xc-1][j-1], ccy[xc-1][j-1], cx[xc-1][j-1], cy[xc-1][j-1], dx[xm-1], dy[k-1]);
    A->a[xm-1].D.el[0][1] = J_gLx_gLy(ccx[xc-1][j-1], ccy[xc-1][j-1], cx[xc-1][j-1], cy[xc-1][j-1], dx[xm-1], dy[k-1]);
    A->a[xm-1].D.el[1][0] = J_gLx_fRy(ccx[xc-1][j],   ccy[xc-1][j],   cx[xc-1][j],   cy[xc-1][j],   dx[xm-1], dy[k]);
    A->a[xm-1].D.el[1][1] = J_gLx_gRy(ccx[xc-1][j],   ccy[xc-1][j],   cx[xc-1][j],   cy[xc-1][j],   dx[xm-1], dy[k]);
    }
else
    {
    A->a[xm-1].B.el[0][0] = J_fLx_fLy(ccx[xc-2][j-1], ccy[xc-2][j-1], cx[xc-2][j-1], cy[xc-2][j-1], dx[xm-1], dy[k-1]);
    A->a[xm-1].B.el[0][1] = J_fLx_gLy(ccx[xc-2][j-1], ccy[xc-2][j-1], cx[xc-2][j-1], cy[xc-2][j-1], dx[xm-1], dy[k-1]);
    A->a[xm-1].B.el[1][0] = J_fLx_fRy(ccx[xc-2][j],   ccy[xc-2][j],   cx[xc-2][j],   cy[xc-2][j],   dx[xm-1], dy[k]);
    A->a[xm-1].B.el[1][1] = J_fLx_gRy(ccx[xc-2][j],   ccy[xc-2][j],   cx[xc-2][j],   cy[xc-2][j],   dx[xm-1], dy[k]);
    
    A->a[xm-1].D.el[0][0] = J_fLx_fLy(ccx[xc-1][j-1], ccy[xc-1][j-1], cx[xc-1][j-1], cy[xc-1][j-1], dx[xm-1], dy[k-1]);
    A->a[xm-1].D.el[0][1] = J_fLx_gLy(ccx[xc-1][j-1], ccy[xc-1][j-1], cx[xc-1][j-1], cy[xc-1][j-1], dx[xm-1], dy[k-1]);
    A->a[xm-1].D.el[1][0] = J_fLx_fRy(ccx[xc-1][j],   ccy[xc-1][j],   cx[xc-1][j],   cy[xc-1][j],   dx[xm-1], dy[k]);
    A->a[xm-1].D.el[1][1] = J_fLx_gRy(ccx[xc-1][j],   ccy[xc-1][j],   cx[xc-1][j],   cy[xc-1][j],   dx[xm-1], dy[k]);
    }

for (i = 0; i <= xm-2; i++)
    {
    A->a[i]  .B    .el[0][0] = J_fLx_fLy(ccx[2*i]  [j-1], ccy[2*i]  [j-1], cx[2*i]  [j-1], cy[2*i]  [j-1], dx[i],   dy[k-1]);
    A->a[i]  .B    .el[0][1] = J_fLx_gLy(ccx[2*i]  [j-1], ccy[2*i]  [j-1], cx[2*i]  [j-1], cy[2*i]  [j-1], dx[i],   dy[k-1]);
    A->a[i]  .B    .el[1][0] = J_fLx_fRy(ccx[2*i]  [j],   ccy[2*i]  [j],   cx[2*i]  [j],   cy[2*i]  [j],   dx[i],   dy[k]);
    A->a[i]  .B    .el[1][1] = J_fLx_gRy(ccx[2*i]  [j],   ccy[2*i]  [j],   cx[2*i]  [j],   cy[2*i]  [j],   dx[i],   dy[k]);
    
    A->a[i]  .D    .el[0][0] = J_fLx_fLy(ccx[2*i+1][j-1], ccy[2*i+1][j-1], cx[2*i+1][j-1], cy[2*i+1][j-1], dx[i],   dy[k-1]);
    A->a[i]  .D    .el[0][1] = J_fLx_gLy(ccx[2*i+1][j-1], ccy[2*i+1][j-1], cx[2*i+1][j-1], cy[2*i+1][j-1], dx[i],   dy[k-1]);
    A->a[i]  .D    .el[1][0] = J_fLx_fRy(ccx[2*i+1][j],   ccy[2*i+1][j],   cx[2*i+1][j],   cy[2*i+1][j],   dx[i],   dy[k]);
    A->a[i]  .D    .el[1][1] = J_fLx_gRy(ccx[2*i+1][j],   ccy[2*i+1][j],   cx[2*i+1][j],   cy[2*i+1][j],   dx[i],   dy[k]);
    
    A->b[i]  .el[0].el[0][0] = J_gLx_fLy(ccx[2*i]  [j-1], ccy[2*i]  [j-1], cx[2*i]  [j-1], cy[2*i]  [j-1], dx[i],   dy[k-1]);
    A->b[i]  .el[0].el[0][1] = J_gLx_gLy(ccx[2*i]  [j-1], ccy[2*i]  [j-1], cx[2*i]  [j-1], cy[2*i]  [j-1], dx[i],   dy[k-1]);
    A->b[i]  .el[0].el[1][0] = J_gLx_fRy(ccx[2*i]  [j],   ccy[2*i]  [j],   cx[2*i]  [j],   cy[2*i]  [j],   dx[i],   dy[k]);
    A->b[i]  .el[0].el[1][1] = J_gLx_gRy(ccx[2*i]  [j],   ccy[2*i]  [j],   cx[2*i]  [j],   cy[2*i]  [j],   dx[i],   dy[k]);
    
    A->b[i]  .el[1].el[0][0] = J_gLx_fLy(ccx[2*i+1][j-1], ccy[2*i+1][j-1], cx[2*i+1][j-1], cy[2*i+1][j-1], dx[i],   dy[k-1]);
    A->b[i]  .el[1].el[0][1] = J_gLx_gLy(ccx[2*i+1][j-1], ccy[2*i+1][j-1], cx[2*i+1][j-1], cy[2*i+1][j-1], dx[i],   dy[k-1]);
    A->b[i]  .el[1].el[1][0] = J_gLx_fRy(ccx[2*i+1][j],   ccy[2*i+1][j],   cx[2*i+1][j],   cy[2*i+1][j],   dx[i],   dy[k]);
    A->b[i]  .el[1].el[1][1] = J_gLx_gRy(ccx[2*i+1][j],   ccy[2*i+1][j],   cx[2*i+1][j],   cy[2*i+1][j],   dx[i],   dy[k]);

    A->c[i]  .el[0].el[0][0] = J_fRx_fLy(ccx[2*i+2][j-1], ccy[2*i+2][j-1], cx[2*i+2][j-1], cy[2*i+2][j-1], dx[i+1], dy[k-1]);
    A->c[i]  .el[0].el[0][1] = J_fRx_gLy(ccx[2*i+2][j-1], ccy[2*i+2][j-1], cx[2*i+2][j-1], cy[2*i+2][j-1], dx[i+1], dy[k-1]);
    A->c[i]  .el[0].el[1][0] = J_fRx_fRy(ccx[2*i+2][j],   ccy[2*i+2][j],   cx[2*i+2][j],   cy[2*i+2][j],   dx[i+1], dy[k]);
    A->c[i]  .el[0].el[1][1] = J_fRx_gRy(ccx[2*i+2][j],   ccy[2*i+2][j],   cx[2*i+2][j],   cy[2*i+2][j],   dx[i+1], dy[k]);
    
    A->c[i]  .el[1].el[0][0] = J_fRx_fLy(ccx[2*i+3][j-1], ccy[2*i+3][j-1], cx[2*i+3][j-1], cy[2*i+3][j-1], dx[i+1], dy[k-1]);
    A->c[i]  .el[1].el[0][1] = J_fRx_gLy(ccx[2*i+3][j-1], ccy[2*i+3][j-1], cx[2*i+3][j-1], cy[2*i+3][j-1], dx[i+1], dy[k-1]);
    A->c[i]  .el[1].el[1][0] = J_fRx_fRy(ccx[2*i+3][j],   ccy[2*i+3][j],   cx[2*i+3][j],   cy[2*i+3][j],   dx[i+1], dy[k]);
    A->c[i]  .el[1].el[1][1] = J_fRx_gRy(ccx[2*i+3][j],   ccy[2*i+3][j],   cx[2*i+3][j],   cy[2*i+3][j],   dx[i+1], dy[k]);
    
    A->a[i+1].A    .el[0][0] = J_gRx_fLy(ccx[2*i+2][j-1], ccy[2*i+2][j-1], cx[2*i+2][j-1], cy[2*i+2][j-1], dx[i+1], dy[k-1]);
    A->a[i+1].A    .el[0][1] = J_gRx_gLy(ccx[2*i+2][j-1], ccy[2*i+2][j-1], cx[2*i+2][j-1], cy[2*i+2][j-1], dx[i+1], dy[k-1]);
    A->a[i+1].A    .el[1][0] = J_gRx_fRy(ccx[2*i+2][j],   ccy[2*i+2][j],   cx[2*i+2][j],   cy[2*i+2][j],   dx[i+1], dy[k]);
    A->a[i+1].A    .el[1][1] = J_gRx_gRy(ccx[2*i+2][j],   ccy[2*i+2][j],   cx[2*i+2][j],   cy[2*i+2][j],   dx[i+1], dy[k]);
    
    A->a[i+1].C    .el[0][0] = J_gRx_fLy(ccx[2*i+3][j-1], ccy[2*i+3][j-1], cx[2*i+3][j-1], cy[2*i+3][j-1], dx[i+1], dy[k-1]);
    A->a[i+1].C    .el[0][1] = J_gRx_gLy(ccx[2*i+3][j-1], ccy[2*i+3][j-1], cx[2*i+3][j-1], cy[2*i+3][j-1], dx[i+1], dy[k-1]);
    A->a[i+1].C    .el[1][0] = J_gRx_fRy(ccx[2*i+3][j],   ccy[2*i+3][j],   cx[2*i+3][j],   cy[2*i+3][j],   dx[i+1], dy[k]);
    A->a[i+1].C    .el[1][1] = J_gRx_gRy(ccx[2*i+3][j],   ccy[2*i+3][j],   cx[2*i+3][j],   cy[2*i+3][j],   dx[i+1], dy[k]);
    }
}  /* end of function load_A_J */


/*************************
 *  FUNCTION:  load_A_L  *
 *************************/
load_A_L(A, j, k, ccx, ccy, cx, cy, dx, dy, t)
struct triblock44 *A;
int j, k;
double **ccx, **ccy, **cx, **cy, *dx, *dy, t;
{
int i;

if (BCTypeW == 1)
    {
    A->a[0].A.el[0][0] = L_gRx_fLy(ccx[0][j-1], ccy[0][j-1], cx[0][j-1], cy[0][j-1], dx[0], dy[k-1] ,t);
    A->a[0].A.el[0][1] = L_gRx_gLy(ccx[0][j-1], ccy[0][j-1], cx[0][j-1], cy[0][j-1], dx[0], dy[k-1] ,t);
    A->a[0].A.el[1][0] = L_gRx_fRy(ccx[0][j],   ccy[0][j],   cx[0][j],   cy[0][j]  , dx[0], dy[k]   ,t);
    A->a[0].A.el[1][1] = L_gRx_gRy(ccx[0][j],   ccy[0][j],   cx[0][j],   cy[0][j]  , dx[0], dy[k]   ,t);
    
    A->a[0].C.el[0][0] = L_gRx_fLy(ccx[1][j-1], ccy[1][j-1], cx[1][j-1], cy[1][j-1], dx[0], dy[k-1] ,t);
    A->a[0].C.el[0][1] = L_gRx_gLy(ccx[1][j-1], ccy[1][j-1], cx[1][j-1], cy[1][j-1], dx[0], dy[k-1] ,t);
    A->a[0].C.el[1][0] = L_gRx_fRy(ccx[1][j],   ccy[1][j],   cx[1][j],   cy[1][j]  , dx[0], dy[k]   ,t);
    A->a[0].C.el[1][1] = L_gRx_gRy(ccx[1][j],   ccy[1][j],   cx[1][j],   cy[1][j]  , dx[0], dy[k]   ,t);
    }
else
    {
    A->a[0].A.el[0][0] = L_fRx_fLy(ccx[0][j-1], ccy[0][j-1], cx[0][j-1], cy[0][j-1], dx[0], dy[k-1] ,t);
    A->a[0].A.el[0][1] = L_fRx_gLy(ccx[0][j-1], ccy[0][j-1], cx[0][j-1], cy[0][j-1], dx[0], dy[k-1] ,t);
    A->a[0].A.el[1][0] = L_fRx_fRy(ccx[0][j],   ccy[0][j],   cx[0][j],   cy[0][j]  , dx[0], dy[k]   ,t);
    A->a[0].A.el[1][1] = L_fRx_gRy(ccx[0][j],   ccy[0][j],   cx[0][j],   cy[0][j]  , dx[0], dy[k]   ,t);
    
    A->a[0].C.el[0][0] = L_fRx_fLy(ccx[1][j-1], ccy[1][j-1], cx[1][j-1], cy[1][j-1], dx[0], dy[k-1] ,t);
    A->a[0].C.el[0][1] = L_fRx_gLy(ccx[1][j-1], ccy[1][j-1], cx[1][j-1], cy[1][j-1], dx[0], dy[k-1] ,t);
    A->a[0].C.el[1][0] = L_fRx_fRy(ccx[1][j],   ccy[1][j],   cx[1][j],   cy[1][j]  , dx[0], dy[k]   ,t);
    A->a[0].C.el[1][1] = L_fRx_gRy(ccx[1][j],   ccy[1][j],   cx[1][j],   cy[1][j]  , dx[0], dy[k]   ,t);
    }

if (BCTypeE == 1)
    {
    A->a[xm-1].B.el[0][0] = L_gLx_fLy(ccx[xc-2][j-1], ccy[xc-2][j-1], cx[xc-2][j-1], cy[xc-2][j-1], dx[xm-1], dy[k-1] ,t);
    A->a[xm-1].B.el[0][1] = L_gLx_gLy(ccx[xc-2][j-1], ccy[xc-2][j-1], cx[xc-2][j-1], cy[xc-2][j-1], dx[xm-1], dy[k-1] ,t);
    A->a[xm-1].B.el[1][0] = L_gLx_fRy(ccx[xc-2][j],   ccy[xc-2][j],   cx[xc-2][j],   cy[xc-2][j],   dx[xm-1], dy[k]   ,t);
    A->a[xm-1].B.el[1][1] = L_gLx_gRy(ccx[xc-2][j],   ccy[xc-2][j],   cx[xc-2][j],   cy[xc-2][j],   dx[xm-1], dy[k]   ,t);
    
    A->a[xm-1].D.el[0][0] = L_gLx_fLy(ccx[xc-1][j-1], ccy[xc-1][j-1], cx[xc-1][j-1], cy[xc-1][j-1], dx[xm-1], dy[k-1] ,t);
    A->a[xm-1].D.el[0][1] = L_gLx_gLy(ccx[xc-1][j-1], ccy[xc-1][j-1], cx[xc-1][j-1], cy[xc-1][j-1], dx[xm-1], dy[k-1] ,t);
    A->a[xm-1].D.el[1][0] = L_gLx_fRy(ccx[xc-1][j],   ccy[xc-1][j],   cx[xc-1][j],   cy[xc-1][j],   dx[xm-1], dy[k]   ,t);
    A->a[xm-1].D.el[1][1] = L_gLx_gRy(ccx[xc-1][j],   ccy[xc-1][j],   cx[xc-1][j],   cy[xc-1][j],   dx[xm-1], dy[k]   ,t);
    }
else
    {
    A->a[xm-1].B.el[0][0] = L_fLx_fLy(ccx[xc-2][j-1], ccy[xc-2][j-1], cx[xc-2][j-1], cy[xc-2][j-1], dx[xm-1], dy[k-1] ,t);
    A->a[xm-1].B.el[0][1] = L_fLx_gLy(ccx[xc-2][j-1], ccy[xc-2][j-1], cx[xc-2][j-1], cy[xc-2][j-1], dx[xm-1], dy[k-1] ,t);
    A->a[xm-1].B.el[1][0] = L_fLx_fRy(ccx[xc-2][j],   ccy[xc-2][j],   cx[xc-2][j],   cy[xc-2][j],   dx[xm-1], dy[k]   ,t);
    A->a[xm-1].B.el[1][1] = L_fLx_gRy(ccx[xc-2][j],   ccy[xc-2][j],   cx[xc-2][j],   cy[xc-2][j],   dx[xm-1], dy[k]   ,t);
    
    A->a[xm-1].D.el[0][0] = L_fLx_fLy(ccx[xc-1][j-1], ccy[xc-1][j-1], cx[xc-1][j-1], cy[xc-1][j-1], dx[xm-1], dy[k-1] ,t);
    A->a[xm-1].D.el[0][1] = L_fLx_gLy(ccx[xc-1][j-1], ccy[xc-1][j-1], cx[xc-1][j-1], cy[xc-1][j-1], dx[xm-1], dy[k-1] ,t);
    A->a[xm-1].D.el[1][0] = L_fLx_fRy(ccx[xc-1][j],   ccy[xc-1][j],   cx[xc-1][j],   cy[xc-1][j],   dx[xm-1], dy[k]   ,t);
    A->a[xm-1].D.el[1][1] = L_fLx_gRy(ccx[xc-1][j],   ccy[xc-1][j],   cx[xc-1][j],   cy[xc-1][j],   dx[xm-1], dy[k]   ,t);
    }

for (i = 0; i <= xm-2; i++)
    {
    A->a[i]  .B    .el[0][0] = L_fLx_fLy(ccx[2*i]  [j-1], ccy[2*i]  [j-1], cx[2*i]  [j-1], cy[2*i]  [j-1], dx[i],   dy[k-1] ,t);
    A->a[i]  .B    .el[0][1] = L_fLx_gLy(ccx[2*i]  [j-1], ccy[2*i]  [j-1], cx[2*i]  [j-1], cy[2*i]  [j-1], dx[i],   dy[k-1] ,t);
    A->a[i]  .B    .el[1][0] = L_fLx_fRy(ccx[2*i]  [j],   ccy[2*i]  [j],   cx[2*i]  [j],   cy[2*i]  [j],   dx[i],   dy[k]   ,t);
    A->a[i]  .B    .el[1][1] = L_fLx_gRy(ccx[2*i]  [j],   ccy[2*i]  [j],   cx[2*i]  [j],   cy[2*i]  [j],   dx[i],   dy[k]   ,t);
    
    A->a[i]  .D    .el[0][0] = L_fLx_fLy(ccx[2*i+1][j-1], ccy[2*i+1][j-1], cx[2*i+1][j-1], cy[2*i+1][j-1], dx[i],   dy[k-1] ,t);
    A->a[i]  .D    .el[0][1] = L_fLx_gLy(ccx[2*i+1][j-1], ccy[2*i+1][j-1], cx[2*i+1][j-1], cy[2*i+1][j-1], dx[i],   dy[k-1] ,t);
    A->a[i]  .D    .el[1][0] = L_fLx_fRy(ccx[2*i+1][j],   ccy[2*i+1][j],   cx[2*i+1][j],   cy[2*i+1][j],   dx[i],   dy[k]   ,t);
    A->a[i]  .D    .el[1][1] = L_fLx_gRy(ccx[2*i+1][j],   ccy[2*i+1][j],   cx[2*i+1][j],   cy[2*i+1][j],   dx[i],   dy[k]   ,t);
    
    A->b[i]  .el[0].el[0][0] = L_gLx_fLy(ccx[2*i]  [j-1], ccy[2*i]  [j-1], cx[2*i]  [j-1], cy[2*i]  [j-1], dx[i],   dy[k-1] ,t);
    A->b[i]  .el[0].el[0][1] = L_gLx_gLy(ccx[2*i]  [j-1], ccy[2*i]  [j-1], cx[2*i]  [j-1], cy[2*i]  [j-1], dx[i],   dy[k-1] ,t);
    A->b[i]  .el[0].el[1][0] = L_gLx_fRy(ccx[2*i]  [j],   ccy[2*i]  [j],   cx[2*i]  [j],   cy[2*i]  [j],   dx[i],   dy[k]   ,t);
    A->b[i]  .el[0].el[1][1] = L_gLx_gRy(ccx[2*i]  [j],   ccy[2*i]  [j],   cx[2*i]  [j],   cy[2*i]  [j],   dx[i],   dy[k]   ,t);
    
    A->b[i]  .el[1].el[0][0] = L_gLx_fLy(ccx[2*i+1][j-1], ccy[2*i+1][j-1], cx[2*i+1][j-1], cy[2*i+1][j-1], dx[i],   dy[k-1] ,t);
    A->b[i]  .el[1].el[0][1] = L_gLx_gLy(ccx[2*i+1][j-1], ccy[2*i+1][j-1], cx[2*i+1][j-1], cy[2*i+1][j-1], dx[i],   dy[k-1] ,t);
    A->b[i]  .el[1].el[1][0] = L_gLx_fRy(ccx[2*i+1][j],   ccy[2*i+1][j],   cx[2*i+1][j],   cy[2*i+1][j],   dx[i],   dy[k]   ,t);
    A->b[i]  .el[1].el[1][1] = L_gLx_gRy(ccx[2*i+1][j],   ccy[2*i+1][j],   cx[2*i+1][j],   cy[2*i+1][j],   dx[i],   dy[k]   ,t);

    A->c[i]  .el[0].el[0][0] = L_fRx_fLy(ccx[2*i+2][j-1], ccy[2*i+2][j-1], cx[2*i+2][j-1], cy[2*i+2][j-1], dx[i+1], dy[k-1] ,t);
    A->c[i]  .el[0].el[0][1] = L_fRx_gLy(ccx[2*i+2][j-1], ccy[2*i+2][j-1], cx[2*i+2][j-1], cy[2*i+2][j-1], dx[i+1], dy[k-1] ,t);
    A->c[i]  .el[0].el[1][0] = L_fRx_fRy(ccx[2*i+2][j],   ccy[2*i+2][j],   cx[2*i+2][j],   cy[2*i+2][j],   dx[i+1], dy[k]   ,t);
    A->c[i]  .el[0].el[1][1] = L_fRx_gRy(ccx[2*i+2][j],   ccy[2*i+2][j],   cx[2*i+2][j],   cy[2*i+2][j],   dx[i+1], dy[k]   ,t);
    
    A->c[i]  .el[1].el[0][0] = L_fRx_fLy(ccx[2*i+3][j-1], ccy[2*i+3][j-1], cx[2*i+3][j-1], cy[2*i+3][j-1], dx[i+1], dy[k-1] ,t);
    A->c[i]  .el[1].el[0][1] = L_fRx_gLy(ccx[2*i+3][j-1], ccy[2*i+3][j-1], cx[2*i+3][j-1], cy[2*i+3][j-1], dx[i+1], dy[k-1] ,t);
    A->c[i]  .el[1].el[1][0] = L_fRx_fRy(ccx[2*i+3][j],   ccy[2*i+3][j],   cx[2*i+3][j],   cy[2*i+3][j],   dx[i+1], dy[k]   ,t);
    A->c[i]  .el[1].el[1][1] = L_fRx_gRy(ccx[2*i+3][j],   ccy[2*i+3][j],   cx[2*i+3][j],   cy[2*i+3][j],   dx[i+1], dy[k]   ,t);
    
    A->a[i+1].A    .el[0][0] = L_gRx_fLy(ccx[2*i+2][j-1], ccy[2*i+2][j-1], cx[2*i+2][j-1], cy[2*i+2][j-1], dx[i+1], dy[k-1] ,t);
    A->a[i+1].A    .el[0][1] = L_gRx_gLy(ccx[2*i+2][j-1], ccy[2*i+2][j-1], cx[2*i+2][j-1], cy[2*i+2][j-1], dx[i+1], dy[k-1] ,t);
    A->a[i+1].A    .el[1][0] = L_gRx_fRy(ccx[2*i+2][j],   ccy[2*i+2][j],   cx[2*i+2][j],   cy[2*i+2][j],   dx[i+1], dy[k]   ,t);
    A->a[i+1].A    .el[1][1] = L_gRx_gRy(ccx[2*i+2][j],   ccy[2*i+2][j],   cx[2*i+2][j],   cy[2*i+2][j],   dx[i+1], dy[k]   ,t);
    
    A->a[i+1].C    .el[0][0] = L_gRx_fLy(ccx[2*i+3][j-1], ccy[2*i+3][j-1], cx[2*i+3][j-1], cy[2*i+3][j-1], dx[i+1], dy[k-1] ,t);
    A->a[i+1].C    .el[0][1] = L_gRx_gLy(ccx[2*i+3][j-1], ccy[2*i+3][j-1], cx[2*i+3][j-1], cy[2*i+3][j-1], dx[i+1], dy[k-1] ,t);
    A->a[i+1].C    .el[1][0] = L_gRx_fRy(ccx[2*i+3][j],   ccy[2*i+3][j],   cx[2*i+3][j],   cy[2*i+3][j],   dx[i+1], dy[k]   ,t);
    A->a[i+1].C    .el[1][1] = L_gRx_gRy(ccx[2*i+3][j],   ccy[2*i+3][j],   cx[2*i+3][j],   cy[2*i+3][j],   dx[i+1], dy[k]   ,t);
    }
}  /* end of function load_A_L */


/********************************************************
 *  FUNCTION:  mult_collvec_collmat -- computes R = MU  *
 *  INPUTS: M = known collmat                           *
 *          U = known vector                            *
 *          R = result of M*U                           *
 ********************************************************/
void mult_collvec_collmat (R, M, U)
struct collvec *R, *U;
struct collmat *M;
{
struct kvec2 tempk1, tempk2;
struct kvec4 k41, k42, k43;
int i, im1;

/********************************
 *  compute RF = AF*UF + BF*U0  *
 ********************************/
mult_M24(&tempk1, &M->BF, &U->M[0]);
mult_tri2xu(&tempk2, &M->AF, &U->F);
addkv2(&R->F, &tempk1, &tempk2);

/****************************************
 *  compute R0 = CF*UF + A0*U0 + B0*U1  *
 ****************************************/
mult_M42(&tempk1, &M->CF, &U->F);
mult_M24(&tempk2, &M->B[0], &U->M[1]);
seam(&k41, &tempk1, &tempk2);
mult_M44(&k42, &M->A[0], &U->M[0]);
addkv4(&R->M[0], &k41, &k42);
    
/************************************************************
 *  compute R(i) = C(i-1)*U(i-1) + A(i)*U(i) + B(i)*U(i+1)  *
 ************************************************************/
for (i = 1; i <= ym-3; i++)
    {
    im1 = i-1;
    mult_M24(&tempk1, &M->C[im1], &U->M[im1]);
    mult_M24(&tempk2, &M->B[i], &U->M[i+1]);
    seam(&k41, &tempk1, &tempk2);
    mult_M44(&k42, &M->A[i], &U->M[i]);
    addkv4(&R->M[i], &k41, &k42);
    }

/*********************************************************************
 *  compute R(ym-2) = C(ym-3)*U(ym-3) + A(ym-2)*U(ym-2) + B(L)*U(L)  *
 *********************************************************************/
i = ym-2;
im1 = i-1;
mult_M24(&tempk1, &M->C[im1], &U->M[im1]);
mult_M42(&tempk2, &M->BL, &U->L);
seam(&k41, &tempk1, &tempk2);
mult_M44(&k42, &M->A[i], &U->M[i]);
addkv4(&R->M[i], &k41, &k42);
    
/*********************************************
 *  compute R(L) = C(L)*U(ym-2) + A(L)*U(L)  *
 *********************************************/
mult_M24(&tempk1, &M->CL, &U->M[ym-2]);
mult_tri2xu(&tempk2, &M->AL, &U->L);
addkv2(&R->L, &tempk1, &tempk2);
}  /* end of function mult_collvec_collmat */


/**********************************************
 *  FUNCTION:  mult_collvec_collmat_trick --  *
 *      computes V = AY   using C             *               
 *  INPUTS: A = known collmat                 *
 *          Y = known vector                  *
 *          V = result of A*Y                 *
 **********************************************/
void mult_collvec_collmat_trick (V, A, Y, C)
struct collvec *V, *Y, *C;
struct collmat *A;
{
struct kvec2 tempk1, tempk2;
struct kvec4 k41, k42, k43;
int i, im1, j, k, l;

/*********
 *  RED  *
 *********/
mult_M24(&tempk1, &A->BF, &Y->M[0]);
addkv2(&V->F, &tempk1, &C->F);

for (i = 1; i <= ym-3; i+=2)
    {
    im1 = i-1;
    mult_M24(&tempk1, &A->C[im1], &Y->M[im1]);
    mult_M24(&tempk2, &A->B[i], &Y->M[i+1]);
    seam(&k41, &tempk1, &tempk2);
    mult_M44(&k42, &A->A[i], &Y->M[i]);
    addkv4(&V->M[i], &k41, &k42);
    }

mult_M24(&tempk1, &A->CL, &Y->M[ym-2]);
addkv2(&V->L, &tempk1, &C->L);

/***********
 *  BLACK  *
 ***********/
for (i = 0; i <= ym-2; i+=2)
    for (j = 0; j < xm; j++)
        for (k = 0; k < 2; k++)
	    for (l = 0; l < 2; l++)
	        V->M[i].el[j].el[k].el[l] = C->M[i].el[j].el[k].el[l];

}  /* end of function mult_collvec_collmat_trick */


/****************************************************
 *  FUNCTION: sub_collvec -- subtract two collvecs  *
 *  INPUTS:  r, a, b   (r = a-b)                    *
 ****************************************************/
void sub_collvec (R, A, B)
struct collvec *R, *A, *B;
{
int i;

subkv2(&R->F, &A->F, &B->F);
for (i = 0; i < ym-1; i++)
    subkv4(&R->M[i], &A->M[i], &B->M[i]);
subkv2(&R->L, &A->L, &B->L);
}  /* end of function sub_collvec */


/****************************************************************
 *  FUNCTION:  copy_collvec -- copies one collvec into another  *
 ****************************************************************/
void copy_collvec (A, B)
struct collvec *A, *B;
{
int i, j, k, l;

for (j = 0; j < xm; j++)
    for (k = 0; k < 2; k++)
        {
	A->F.el[j].el[k] = B->F.el[j].el[k];
	A->L.el[j].el[k] = B->L.el[j].el[k];
        }

for (i = 0; i <= ym-2; i++)
    for (j = 0; j < xm; j++)
        for (k = 0; k < 2; k++)
	    for (l = 0; l < 2; l++)
	        A->M[i].el[j].el[k].el[l] = B->M[i].el[j].el[k].el[l];
}  /* end of function copy_collvec */


/************************************************************************
 *  FUNCTION:  dot_collvec --determine the dot product of two collvecs  *
 *  INPUTS: A, B                                                        *
 *  OUTPUT: A (dot) B                                                   *
 ************************************************************************/
double dot_collvec (A, B)
struct collvec *A, *B;
{
int i, j, k, l;
double sum = 0;

for (i = 0; i < xm; i++)
    for (j = 0; j < 2; j++)
        {
	sum += A->F.el[i].el[j] * B->F.el[i].el[j];
	sum += A->L.el[i].el[j] * B->L.el[i].el[j];
        }

for (i = 0; i <= ym-2; i++)
    for (j = 0; j < xm; j++)
        for (k = 0; k < 2; k++)
	    for (l = 0; l < 2; l++)
	        sum += A->M[i].el[j].el[k].el[l] * B->M[i].el[j].el[k].el[l];

return(sum);
}  /* end of function dot_collvec */


/*******************************************************
 *  FUNCTION:  saxpy_collvec -- compute r=ax+y, where  *
 *             x and y are collvecs                    *
 *             a = scalar                              *
 *******************************************************/
void saxpy_collvec (r, a, x, y)
struct collvec *r, *x, *y;
double a;
{
struct collvec tmp;

scmult_collvec(&tmp, x, a);
add_collvec(r, &tmp, y);
}  /* end of function saxpy_collvec *


/*****************************************************
 *  FUNCTION: psolve -- solves matrix equation Ax=b  *
 *         A = D-L (i.e., A = Gauss-Seidel matrix)   *
 *****************************************************/
void psolve (x, A, b)
struct collvec *x, *b;
struct collmat *A;
{
struct collvec tmp;
struct kvec2 tempk1, tempk2;
struct kvec4 k41;
int i, im1;

/*********
 *  RED  *
 *********/
tridiag2(&x->F, &A->AF, &b->F);

for (i = 1; i <= ym-3; i += 2)
    tridiag4(&x->M[i], &A->A[i], &b->M[i]);

tridiag2(&x->L, &A->AL, &b->L);

/***********
 *  BLACK  *
 ***********/
for (i = 0; i <= ym-2; i += 2)
    {
    if (i == 0)
        {
	mult_M42(&tempk1, &A->CF, &x->F);
	mult_M24(&tempk2, &A->B[0], &x->M[1]);
	seam(&k41, &tempk1, &tempk2);
	subkv4(&tmp.M[0], &b->M[0], &k41);
        }
    else if (i == ym-2)
        {
	im1 = i-1;
	mult_M24(&tempk1, &A->C[im1], &x->M[im1]);
	mult_M42(&tempk2, &A->BL, &x->L);
	seam(&k41, &tempk1, &tempk2);
	subkv4(&tmp.M[i], &b->M[i], &k41);
        }
    else
        {
	im1 = i-1;
	mult_M24(&tempk1, &A->C[im1], &x->M[im1]);
	mult_M24(&tempk2, &A->B[i], &x->M[i+1]);
	seam(&k41, &tempk1, &tempk2);
	subkv4(&tmp.M[i], &b->M[i], &k41);
        }

    tridiag4(&x->M[i], &A->A[i], &tmp.M[i]);
    }
}  /* end of function psolve */


/********************************************************************
 *  FUNCTION:  norm_collvec -- find the infinity norm of a collvec  *
 *  INPUT:  A = collvec                                             *
 *  OUTPUT: nm = norm                                               *
 ********************************************************************/
double norm_collvec (A)
struct collvec *A;
{
int i, j, k, l;
double nm = 0;

for (i = 0; i < xm; i++)
    for (j = 0; j < 2; j++)
        {
	if (fabs(A->F.el[i].el[j]) > nm)
	    nm = fabs(A->F.el[i].el[j]);
	if (fabs(A->L.el[i].el[j]) > nm)
	    nm = fabs(A->L.el[i].el[j]);
        }

for (i = 0; i < ym-1; i++)
    for (j = 0; j < xm; j++)
        for (k = 0; k < 2; k++)
	    for (l = 0; l < 2; l++)
	        if (fabs(A->M[i].el[j].el[k].el[l]) > nm)
		    nm = fabs(A->M[i].el[j].el[k].el[l]);
return(nm);
}  /* end of function norm_collvec */