calc_interp_weights.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_INTERP_WEIGHTS                                     !
!                                                                      !
!  Purpose: Calculate weights used for interpolation.                  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_interp_weights

      use particle_filter, only: des_interp_scheme_enum
      use particle_filter, only: des_interp_dpvm
      use particle_filter, only: des_interp_gauss
      use particle_filter, only: des_interp_lhat

      SELECT CASE(des_interp_scheme_enum)
      CASE(des_interp_dpvm);  CALL calc_interp_weights1
      CASE(des_interp_gauss); CALL calc_interp_weights1
      CASE(des_interp_lhat);  CALL calc_interp_weights2
      END SELECT

      RETURN
      END SUBROUTINE calc_interp_weights

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_INTERP_WEIGHTS1                                    !
!                                                                      !
!  Purpose: Calculate weights used for interpolation.                  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_interp_weights1

      use discretelement, only: max_pip
      use discretelement, only: pijk
      use discretelement, only: des_pos_new
      use discretelement, only: xe, yn, zt
      use particle_filter, only: filter_cell
      use particle_filter, only: filter_weight
      use geometry, only: do_k

      use param1, only: zero, one
! compar is needed by functions.inc
      use compar

      IMPLICIT NONE

      INTEGER :: L, IDX
      INTEGER :: IJK, IJKt

      INTEGER :: I, J, K
      INTEGER :: IC, JC, KC
      INTEGER :: Km, Kp

      DOUBLE PRECISION :: WEIGHT
      DOUBLE PRECISION :: WEIGHT_I(-1:1)
      DOUBLE PRECISION :: WEIGHT_J(-1:1)
      DOUBLE PRECISION :: WEIGHT_K(-1:1)

!$omp parallel default(none) private(L, WEIGHT_I, WEIGHT_J, WEIGHT_K,  &
!$omp    KM, KP, IDX, IJK, I, J, K, WEIGHT, IJKT)                      &
!$omp shared(MAX_PIP, PIJK, DES_POS_NEW, XE, YN, ZT, DO_K,        &
!$omp    FILTER_CELL, FILTER_WEIGHT)
!$omp do
      DO l = 1, max_pip

         IF(.NOT.is_normal(l)) cycle

         i = pijk(l,1)
         j = pijk(l,2)
         k = pijk(l,3)

! Tentative weights for I indices to the West and East.
         weight_i(-1) = calc_filter_weights(des_pos_new(l,1), xe(i-1))
         weight_i( 1) = calc_filter_weights(xe(i), des_pos_new(l,1))
         weight_i( 0) = one - weight_i(-1) - weight_i(1)

! Tentative weights for J indices to the South and North.
         weight_j(-1) = calc_filter_weights(des_pos_new(l,2), yn(j-1))
         weight_j( 1) = calc_filter_weights(yn(j), des_pos_new(l,2))
         weight_j( 0) = one - weight_j(-1) - weight_j(1)

! Tentative weights for K indices to the Top and Bottom.
         IF(do_k) THEN
            km=-1;  kp=1
            weight_k(-1) = calc_filter_weights(des_pos_new(l,3),zt(k-1))
            weight_k( 1) = calc_filter_weights(zt(k), des_pos_new(l,3))
            weight_k( 0) = one - weight_k(-1) - weight_k(1)
         ELSE
            km= 0; kp=0
            weight_k( 0) = one
         ENDIF

! Set the default fluid cell index and loop counter.
         ijk = pijk(l,4)
         idx=0

! Calculate weights for ghost particles. Only store weights that the
! current process owns.
         DO kc=km,kp
         DO ic=-1,+1
         DO jc=-1,+1
            idx=idx+1
            weight = weight_i(ic)*weight_j(jc)*weight_k(kc)
            IF(is_on_mype_plus2layers(i+ic,j+jc,k+kc)) THEN
               ijkt = funijk(i+ic,j+jc,k+kc)
               IF(fluid_at(ijkt) .OR. cyclic_at(ijkt)) THEN
                  filter_cell(idx,l) = ijkt
                  filter_weight(idx,l) = weight
               ELSE
                  filter_cell(idx,l) = ijk
                  filter_weight(idx,l) = weight
               ENDIF
            ELSE
               filter_cell(idx,l) = ijk
               filter_weight(idx,l) = zero
            ENDIF
         ENDDO
         ENDDO
         ENDDO

      ENDDO
!$omp end do
!$omp end parallel

      CONTAINS

!``````````````````````````````````````````````````````````````````````!
!                                                                      !
!``````````````````````````````````````````````````````````````````````!
      DOUBLE PRECISION FUNCTION calc_filter_weights(POS1, POS2)

      use particle_filter, only: filter_width_interp
      use particle_filter, only: oofilter_vol
      use particle_filter, only: filter_width_interpx3

      DOUBLE PRECISION, INTENT(IN) :: POS1, POS2
      DOUBLE PRECISION :: OVERLAP, HEIGHT

      overlap = pos1 - pos2
      IF(overlap < filter_width_interp) THEN
         height = filter_width_interp - overlap
         calc_filter_weights = height**2 * &
            (filter_width_interpx3 - height)*oofilter_vol
      ELSE
         calc_filter_weights = zero
      ENDIF

      END FUNCTION calc_filter_weights

      include 'functions.inc'

      END SUBROUTINE calc_interp_weights1


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_INTERP_WEIGHTS2                                    !
!                                                                      !
!  Purpose: Calculate weights used for interpolation.                  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_interp_weights2

      use discretelement, only: max_pip
      use discretelement, only: pijk
      use discretelement, only: des_pos_new
      use discretelement, only: xe, yn, zt
      use particle_filter, only: filter_cell
      use particle_filter, only: filter_weight
      use geometry, only: dx, dy, dz
      use geometry, only: odx_e, ody_n, odz_t, do_k

      use param1, only: zero, half, one

      use compar, only: funijk_map_c


      IMPLICIT NONE

      INTEGER :: L, IDX
      INTEGER :: IJK, IJKt

      INTEGER :: I, J, K
      INTEGER :: IC, JC, KC
      INTEGER :: Km, Kp

      DOUBLE PRECISION :: XC, YC, ZC
      DOUBLE PRECISION :: WEIGHT
      DOUBLE PRECISION :: WEIGHT_I(-1:1)
      DOUBLE PRECISION :: WEIGHT_J(-1:1)
      DOUBLE PRECISION :: WEIGHT_K(-1:1)


      DO l = 1, max_pip

         IF(.NOT.is_normal(l)) cycle

         i = pijk(l,1)
         j = pijk(l,2)
         k = pijk(l,3)

! Tentative weights for I indices to the West and East.
         xc = xe(i-1) + half*dx(i)
         IF(des_pos_new(l,1) < xe(i-1)+half*dx(i)) THEN
            weight_i(-1) = (xc - des_pos_new(l,1))*odx_e(i-1)
            weight_i( 0) = one - weight_i(-1)
            weight_i( 1) = zero
         ELSE
            weight_i( 1) = (des_pos_new(l,1) - xc)*odx_e(i)
            weight_i( 0) = one - weight_i(1)
            weight_i(-1) = zero
         ENDIF

! Tentative weights for J indices to the South and North.
         yc = yn(j-1) + half*dy(j)
         IF(des_pos_new(l,2) < yn(j-1)+half*dy(j)) THEN
            weight_j(-1) = (yc - des_pos_new(l,2))*ody_n(j-1)
            weight_j( 0) = one - weight_j(-1)
            weight_j( 1) = zero
         ELSE
            weight_j( 1) = (des_pos_new(l,2) - yc)*ody_n(j)
            weight_j( 0) = one - weight_j(1)
            weight_j(-1) = zero
         ENDIF

! Tentative weights for K indices to the Top and Bottom.
         IF(do_k) THEN
            km=-1;  kp=1
            zc = zt(k-1) + half*dz(k)
            IF(des_pos_new(l,3) < zt(k-1)+half*dz(k)) THEN
               weight_k(-1) = (zc - des_pos_new(l,3))*odz_t(k-1)
               weight_k( 0) = one - weight_k(-1)
               weight_k( 1) = zero
            ELSE
               weight_k( 1) = (des_pos_new(l,3) - zc)*odz_t(k)
               weight_k( 0) = one - weight_k(1)
               weight_k(-1) = zero
            ENDIF
         ELSE
            km= 0; kp=0
            weight_k( 0) = one
         ENDIF

! Set the default fluid cell index and loop counter.
         ijk = pijk(l,4)
         idx=0

! Calculate weights for ghost particles. Only store weights that the
! current process owns.
         DO kc=km,kp
         DO ic=-1,+1
         DO jc=-1,+1
            idx=idx+1
            weight = weight_i(ic)*weight_j(jc)*weight_k(kc)
            IF(is_on_mype_plus2layers(i+ic,j+jc,k+kc)) THEN
               ijkt = funijk_map_c(i+ic,j+jc,k+kc)
               IF(fluid_at(ijkt)) THEN
                  filter_cell(idx,l) = ijkt
                  filter_weight(idx,l) = weight
               ELSE
                  filter_cell(idx,l) = ijk
                  filter_weight(idx,l) = weight
               ENDIF
            ELSE
               filter_cell(idx,l) = ijk
               filter_weight(idx,l) = zero
            ENDIF
         ENDDO
         ENDDO
         ENDDO

      ENDDO

      CONTAINS

      include 'functions.inc'

      END SUBROUTINE calc_interp_weights2