apply_gradient_bc_face Subroutine

    subroutine apply_gradient_bc_face(qp, temp, cells, faces, dims, domain, bc_id, fixed_temp)
      !< Call same subroutine for all the face
      !< Apply/set value of all gradient in the ghost cells
      !< gradqp_G = (qp_I - qp_G)*Area_W*unit_normal_G/(volume_G)
      !< volume_G = volume_I
      !-----------------------------------------------------------
      implicit none
      type(extent), intent(in) :: dims
      !< Extent of the domain:imx,jmx,kmx
      type(singlesub), intent(in) :: domain
      !< flags for direction 
      real(wp), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2, 1:dims%n_var), intent(in) :: qp
      !< Input variable of which graident is required
      real(wp), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2), intent(in) :: Temp
      !< Intput Temperature variable
      type(celltype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: cells
      !< Input cell quantities: volume
      type(facetype), dimension(-2:dims%imx+3,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: faces
      !< Input varaible which stores any(I,J,K) faces' area and unit normal
      integer, intent(in) :: bc_id
      real(wp), intent(in) :: fixed_temp
      real(wp), dimension(n_grad) :: qp_I
      real(wp), dimension(n_grad) :: qp_G
      real(wp)    :: T_I
      real(wp)    :: T_G 
      real(wp)    :: c_x
      real(wp)    :: c_y
      real(wp)    :: c_z
      integer :: i,j,k,l
      real(wp)    :: nx
      real(wp)    :: ny
      real(wp)    :: nz
      real(wp)    :: dot
      integer  :: il, jl, kl
      integer  :: iu, ju, ku

      il = domain%il
      jl = domain%jl
      kl = domain%kl
      iu = domain%iu
      ju = domain%ju
      ku = domain%ku

      do k = domain%kmin, domain%kmax
        do j = domain%jmin, domain%jmax
          do i = domain%imin, domain%imax
            nx   = faces(i,j,k)%nx
            ny   = faces(i,j,k)%ny
            nz   = faces(i,j,k)%nz
            c_x  = faces(i,j,k)%A*nx/cells(i-iu,j-ju,k-ku)%volume
            c_y  = faces(i,j,k)%A*ny/cells(i-iu,j-ju,k-ku)%volume
            c_z  = faces(i,j,k)%A*nz/cells(i-iu,j-ju,k-ku)%volume
            T_I  = Temp(i-iu,j-ju,k-ku)
            T_G  = Temp(i-il,j-jl,k-kl)
            qp_I = qp(i-iu,j-ju,k-ku,2:dims%n_var)
            qp_G = qp(i-il,j-jl,k-kl,2:dims%n_var)

            ! normal component of gradient
            gradqp_x(i-il,j-jl,k-kl,:) = domain%sig*(qp_I - qp_G)*c_x 
            gradqp_y(i-il,j-jl,k-kl,:) = domain%sig*(qp_I - qp_G)*c_y
            gradqp_z(i-il,j-jl,k-kl,:) = domain%sig*(qp_I - qp_G)*c_z
            gradqp_x(i-il,j-jl,k-kl,4) = domain%sig*( T_I -  T_G)*c_x
            gradqp_y(i-il,j-jl,k-kl,4) = domain%sig*( T_I -  T_G)*c_y
            gradqp_z(i-il,j-jl,k-kl,4) = domain%sig*( T_I -  T_G)*c_z
            if(bc_id==-5 .and. (fixed_temp<1. .and. fixed_temp>=0.))then
              gradqp_x(i-il,j-jl,k-kl,4) = -gradqp_x(i-iu,j-ju,k-ku,4)
              gradqp_y(i-il,j-jl,k-kl,4) = -gradqp_y(i-iu,j-ju,k-ku,4)
              gradqp_z(i-il,j-jl,k-kl,4) = -gradqp_z(i-iu,j-ju,k-ku,4)

            end if
            !parallel component of gradient
            do l=1,n_grad
              dot = (gradqp_x(i-iu,j-ju,k-ku,l)*nx) + (gradqp_y(i-iu,j-ju,k-ku,l)*ny) + (gradqp_z(i-iu,j-ju,k-ku,l)*nz)
              gradqp_x(i-il,j-jl,k-kl,l) = gradqp_x(i-il,j-jl,k-kl,l) + (gradqp_x(i-iu,j-ju,k-ku,l) - dot*nx)
              gradqp_y(i-il,j-jl,k-kl,l) = gradqp_y(i-il,j-jl,k-kl,l) + (gradqp_y(i-iu,j-ju,k-ku,l) - dot*ny)
              gradqp_z(i-il,j-jl,k-kl,l) = gradqp_z(i-il,j-jl,k-kl,l) + (gradqp_z(i-iu,j-ju,k-ku,l) - dot*nz)
            end do
          end do
        end do
      end do

    end subroutine apply_gradient_bc_face