dif_w_is.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DIF_W_IS                                                C
!  Purpose: Remove diffusive fluxes across internal surfaces.          C
!  (Make user defined internal surfaces non-conducting)                C
!                                                                      C
!  Author: M. Syamlal                                 Date: 30-APR-97  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE dif_w_is(DIF, A_M, M)

! Modules
!---------------------------------------------------------------------//
      USE param
      USE geometry, only: odx_e, ody_n, axz_w, ayz_w

      USE is, only: is_defined, is_plane
      USE is, only: is_i_w, is_i_e, is_j_s, is_j_n, is_k_t, is_k_b

      USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h

      USE functions, only: funijk, ip_of, jp_of
      USE functions, only: east_of, north_of, top_of
      USE functions, only: is_on_mype_plus2layers

      USE compar, only: dead_cell_at
      IMPLICIT NONE
! Dummy arguments
!---------------------------------------------------------------------//
! Gamma -- diffusion coefficient
      DOUBLE PRECISION, INTENT(IN) :: Dif(dimension_3)
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Solids phase
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
! Diffusion parameter
      DOUBLE PRECISION :: D_f
! Internal surface
      INTEGER :: L
! Indices
      INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK
      INTEGER :: IJKE, IJKN, IJKT, IJPK, IPJK, IJKTN, IJKTE
!---------------------------------------------------------------------//


      DO l = 1, dimension_is
         IF (is_defined(l)) THEN
            i1 = is_i_w(l)
            i2 = is_i_e(l)
            j1 = is_j_s(l)
            j2 = is_j_n(l)
            k1 = is_k_b(l)
            k2 = is_k_t(l)

            IF (is_plane(l) == 'N') THEN
               DO k = k1, k2
               DO j = j1, j2
               DO i = i1, i2
                  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
                  IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
                  ijk = funijk(i,j,k)
                  ijkt = top_of(ijk)
                  ijkn = north_of(ijk)
                  ijktn = top_of(ijkn)
                  ijpk = jp_of(ijk)

                  d_f = avg_z_h(avg_y_h(dif(ijk),dif(ijkn),j),&
                                avg_y_h(dif(ijkt),dif(ijktn),j),k)*&
                        ody_n(j)*axz_w(ijk)

                  a_m(ijk,north,m) = a_m(ijk,north,m) - d_f
                  a_m(ijpk,south,m) = a_m(ijpk,south,m) - d_f
               ENDDO
               ENDDO
               ENDDO

            ELSEIF (is_plane(l) == 'E') THEN
               DO k = k1, k2
               DO j = j1, j2
               DO i = i1, i2
                  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
                  IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
                  ijk = funijk(i,j,k)
                  ijke = east_of(ijk)
                  ijkt = top_of(ijk)
                  ijkte = east_of(ijkt)
                  ipjk = ip_of(ijk)

                  d_f = avg_z_h(avg_x_h(dif(ijk),dif(ijke),i),&
                                avg_x_h(dif(ijkt),dif(ijkte),i),k)*&
                        odx_e(i)*ayz_w(ijk)

                  a_m(ijk,east,m) = a_m(ijk,east,m) - d_f
                  a_m(ipjk,west,m) = a_m(ipjk,west,m) - d_f
               ENDDO
               ENDDO
               ENDDO
            ENDIF

         ENDIF   ! end if is_defined
      ENDDO   ! end do dimension_is
      RETURN
      END SUBROUTINE dif_w_is