/* 
 * ==============================bc====================================
 *    			Boundary Conditions
 *    ** Be sure to modify this routine to be consistent w/your code **
 */

#include "program6_BC.h" 

void bc(float t1[][NYDIM  ][NZDIM  ],float t2[][NYDIM  ][NZDIM  ],
	float u1[][NYDIM  ][NZDIM  ],float u2[][NYDIM  ][NZDIM  ],
	float v1[][NYDIM+1][NZDIM  ],float v2[][NYDIM+1][NZDIM  ],
	float w1[][NYDIM  ][NZDIM+1],float w2[][NYDIM  ][NZDIM+1],
	float p1[][NYDIM  ][NZDIM  ],float p2[][NYDIM  ][NZDIM  ])
{
    int i,j,k;    
    /* X  faces */
    #pragma omp parallel for shared(u1,u2,p1,p2,t1,t2) private(j,k)
    for(j=J1;j<=J2;j++) {
        for (k=K1;k<=K2;k++) {         	    
            u1[I1][j][k]=-u1[I1+1][j][k];
            u2[I1][j][k]=-u2[I1+1][j][k];
            u1[I2+1][j][k]=-u1[I2][j][k];
            u2[I2+1][j][k]=-u2[I2][j][k];
            p1[I1-1][j][k]=p1[I1+1][j][k];
            p2[I1-1][j][k]=p2[I1+1][j][k];
            p1[I2+1][j][k]=p1[I2-1][j][k];
            p2[I2+1][j][k]=p2[I2-1][j][k];
            t1[I1-1][j][k]=t1[I1+1][j][k];
            t2[I1-1][j][k]=t2[I1+1][j][k];
            t1[I2+1][j][k]=t1[I2-1][j][k];
            t2[I2+1][j][k]=t2[I2-1][j][k];
	    /* Need 2nd ghost point for t1 - Brian */
        }
    }
    #pragma omp parallel for shared(v1,v2) private(j,k)
    for(j=J1;j<=J2+1;j++) {
        for (k=K1;k<=K2;k++) {
            v1[I1-1][j][k]=v1[I1+1][j][k];
            v2[I1-1][j][k]=v2[I1+1][j][k];
            v1[I2+1][j][k]=v1[I2-1][j][k];
            v2[I2+1][j][k]=v2[I2-1][j][k];
        }
    }
    #pragma omp parallel for shared(w1,w2) private(j,k)
    for(j=J1;j<=J2;j++) {
        for (k=K1+1;k<=K2;k++) {
            w1[I1-1][j][k]=w1[I1+1][j][k];
            w2[I1-1][j][k]=w2[I1+1][j][k];
            w1[I2+1][j][k]=w1[I2-1][j][k];
            w2[I2+1][j][k]=w2[I2-1][j][k];
        }
    }
	
    /* Y Faces */
    #pragma omp parallel for shared(v1,v2,p1,p2,t1,t2) private(i,k)
    for(i=I1;i<=I2;i++) {
        for(k=K1;k<=K2;k++) {
	    v1[i][J2+1][k]=v1[i][J1][k];
	    v2[i][J2+1][k]=v2[i][J1][k];
	    v1[i][J1-1][k]=v1[i][J2][k];
	    v2[i][J1-1][k]=v2[i][J2][k];
	    t1[i][J2+1][k]=t1[i][J1][k];
	    t2[i][J2+1][k]=t2[i][J1][k];
	    t1[i][J1-1][k]=t1[i][J2][k];
	    t2[i][J1-1][k]=t2[i][J2][k];
	    /* Need 2nd ghost point for t1 - Brian */
	    p1[i][J2+1][k]=p1[i][J1][k];
	    p2[i][J2+1][k]=p2[i][J1][k];
	    p1[i][J1-1][k]=p1[i][J2][k];
	    p2[i][J1-1][k]=p2[i][J2][k];
	}
    }
	
    #pragma omp parallel for shared(u1,u2) private(i,k) 
    for(i=I1;i<=I2+1;i++) {
	for(k=K1;k<=K2;k++) {
	    u1[i][J2+1][k]=u1[i][J1][k];
	    u2[i][J2+1][k]=u2[i][J1][k];
	    u1[i][J1-1][k]=u1[i][J2][k];
	    u2[i][J1-1][k]=u2[i][J2][k];
	}
    }
    #pragma omp parallel for shared(w1,w2) private(i,k) 
    for(i=I1;i<=I2;i++) {
	for(k=K1+1;k<=K2;k++) {
	    w1[i][J2+1][k]=w1[i][J1][k];
	    w2[i][J2+1][k]=w2[i][J1][k];
	    w1[i][J1-1][k]=w1[i][J2][k];
	    w2[i][J1-1][k]=w2[i][J2][k];
	}
    }
	
    /* Z Faces */
    #pragma omp parallel for shared(w1,w2,p1,p2,t1,t2) private(i,j) 
    for(i=I1;i<=I2;i++) {
	for(j=J1;j<=J2;j++) {
            t1[i][j][K1-1] = t1[i][j][K1];
            t2[i][j][K1-1] = t2[i][j][K1];
            t1[i][j][K2+1] = t1[i][j][K2];
            t2[i][j][K2+1] = t2[i][j][K2];
	    /* Need 2nd ghost point for t1 - Brian */
            p1[i][j][K1-1] = p1[i][j][K1];
            p2[i][j][K1-1] = p2[i][j][K1];
            p1[i][j][K2+1] = p1[i][j][K2];
            p2[i][j][K2+1] = p2[i][j][K2];
            w1[i][j][K1]  = 0.0;
            w2[i][j][K1]  = 0.0;
            w1[i][j][K2+1]= 0.0;
            w2[i][j][K2+1]= 0.0;
                  
        }
    }
    #pragma omp parallel for shared(u1,u2) private(i,j)         
    for (i=I1;i<=I2+1;i++) {
        for(j=J1;j<=J2;j++) {
            u1[i][j][K1-1]= u1[i][j][K1];
            u2[i][j][K1-1]= u2[i][j][K1];
            u1[i][j][K2+1]= u1[i][j][K2];
            u2[i][j][K2+1]= u2[i][j][K2];
	}
    }
	#pragma omp parallel for shared(v1,v2) private(i,j) 
	for (i=I1;i<=I2;i++) {
	    for(j=J1;j<=J2+1;j++) {
            v1[i][j][K1-1]= v1[i][j][K1];
            v2[i][j][K1-1]= v2[i][j][K1];
            v1[i][j][K2+1]= v1[i][j][K2];
            v2[i][j][K2+1]= v2[i][j][K2];
		}
    }

      return;
}