dif_u_is.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DIF_U_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_u_is(DIF, A_M, M)

! Modules
!---------------------------------------------------------------------//
      USE param
      USE geometry, only: do_k
      USE geometry, only: ody_n, ox_e, odz_t, axz_u, axy_u

      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, jp_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
      USE compar, only: istart2, jstart2, kstart2
      USE compar, only: iend2, jend2, kend2
      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, IJPK, IJKP, IJKNE, IJKTE
! Diffusion parameter
      DOUBLE PRECISION :: D_f
!---------------------------------------------------------------------//

      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)

! Limit I1, I2 and all to local processor first ghost layer
            IF(i1.LE.iend2)   i1 = max(i1, istart2)
            IF(j1.LE.jend2)   j1 = max(j1, jstart2)
            IF(k1.LE.kend2)   k1 = max(k1, kstart2)
            IF(i2.GE.istart2) i2 = min(i2, iend2)
            IF(j2.GE.jstart2) j2 = min(j2, jend2)
            IF(k2.GE.kstart2) k2 = min(k2, kend2)

            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)
                  ijke = east_of(ijk)
                  ijkn = north_of(ijk)
                  ijkne = east_of(ijkn)
                  ijpk = jp_of(ijk)

                  d_f = avg_x_h(avg_y_h(dif(ijk),dif(ijkn),j),&
                                avg_y_h(dif(ijke),dif(ijkne),j),i)*&
                        ody_n(j)*axz_u(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) == '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)
                     ijke = east_of(ijk)
                     ijkt = top_of(ijk)
                     ijkte = east_of(ijkt)
                     ijkp = kp_of(ijk)

                     d_f = avg_x_h(avg_z_h(dif(ijk),dif(ijkt),k),&
                                   avg_z_h(dif(ijke),dif(ijkte),k),i)*&
                           ox_e(i)*odz_t(k)*axy_u(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_u_is