dif_v_is.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DIF_V_IS                                                C
!  Purpose: Remove diffusive fluxes across internal surfaces.          C
!                                                                      C
!  Author: M. Syamlal                                 Date: 30-APR-97  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE dif_v_is(DIF, A_M, M)

! Modules
!---------------------------------------------------------------------//
      USE param
      USE geometry, only: do_k
      USE geometry, only: odx_e, ox, odz_t, axy_v, ayz_v

      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, kp_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
!---------------------------------------------------------------------//
! Internal surface
      INTEGER :: L
! Indices
      INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK
      INTEGER :: IJKE, IJKN, IJKT, IPJK, IJKP, IJKNE, IJKTN
! Diffusion parameter
      DOUBLE PRECISION :: D_f
!---------------------------------------------------------------------//


! Make user defined internal surfaces non-conducting
      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) == '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)
                  ijkn = north_of(ijk)
                  ijkne = east_of(ijkn)
                  ipjk = ip_of(ijk)

                  d_f = avg_y_h(avg_x_h(dif(ijk),dif(ijke),i),&
                                avg_x_h(dif(ijkn),dif(ijkne),i),j)*&
                        odx_e(i)*ayz_v(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

            ELSEIF (is_plane(l) == 'T') THEN
               IF (do_k) 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)
                     ijkn = north_of(ijk)
                     ijkt = top_of(ijk)
                     ijktn = north_of(ijkt)
                     ijkp = kp_of(ijk)

                     d_f = avg_y_h(avg_z_h(dif(ijk),dif(ijkt),k),&
                                   avg_z_h(dif(ijkn),dif(ijktn),k),j)*&
                           ox(i)*odz_t(k)*axy_v(ijk)

                     a_m(ijk,top,m) = a_m(ijk,top,m) - d_f
                     a_m(ijkp,bottom,m) = a_m(ijkp,bottom,m) - d_f
                  ENDDO
                  ENDDO
                  ENDDO
               ENDIF
            ENDIF

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