diff --git a/sources/boundaries.F90 b/sources/boundaries.F90 index 35988b6..96f0201 100644 --- a/sources/boundaries.F90 +++ b/sources/boundaries.F90 @@ -946,7 +946,7 @@ module boundaries ! prolong corner boundaries from blocks at lower levels ! - call boundaries_corner_prolong() + call boundaries_corner_prolong(status) end if ! minlev /= maxlev @@ -4388,7 +4388,7 @@ module boundaries ! !=============================================================================== ! - subroutine boundaries_corner_prolong() + subroutine boundaries_corner_prolong(status) ! import external procedures and variables ! @@ -4410,6 +4410,8 @@ module boundaries ! implicit none + integer, intent(out) :: status + ! local pointers ! type(block_meta), pointer :: pmeta, pneigh @@ -4529,12 +4531,12 @@ module boundaries #if NDIMS == 2 call block_corner_prolong((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & - , pmeta%data%q(1:nv,il:iu,jl:ju, : )) + , pmeta%data%q(1:nv,il:iu,jl:ju, : ), status) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_corner_prolong((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & - , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku)) + , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku), status) #endif /* NDIMS == 3 */ #ifdef MPI @@ -4644,12 +4646,12 @@ module boundaries #if NDIMS == 2 call block_corner_prolong((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & - , sbuf(l,1:nv,1:ng,1:ng, : )) + , sbuf(l,1:nv,1:ng,1:ng, : ), status) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_corner_prolong((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & - , sbuf(l,1:nv,1:ng,1:ng,1:ng)) + , sbuf(l,1:nv,1:ng,1:ng,1:ng), status) #endif /* NDIMS == 3 */ ! associate the pointer with the next block @@ -5996,38 +5998,45 @@ module boundaries ! subroutine BLOCK_CORNER_PROLONG: ! ------------------------------- ! -! Subroutine returns the corner boundary region by prolongating -! the corresponding region from the provided input variable array. +! Subroutine prolongs the region of qn specified by the corner position +! and inserts it in the corresponding corner boundary region of qb. +! +! Remarks: +! +! This subroutine requires the ghost zone regions of the qn already +! updated in order to perform a consistent interpolation. ! ! Arguments: ! -! pos - the corner position; -! qn - the input neighbor variable array; -! qb - the output corner boundary array; +! pos - the corner position; +! qn - the array to prolong from; +! qb - the array of the corner boundary region; +! status - the call status; ! !=============================================================================== ! - subroutine block_corner_prolong(pos, qn, qb) + subroutine block_corner_prolong(pos, qn, qb, status) - use coordinates , only : corners_dp + use coordinates , only : nb, ne, nh => nghosts_half use equations , only : nv, positive use helpers , only : print_message use interpolations, only : limiter_prol implicit none - integer , dimension(3) , intent(in) :: pos - real(kind=8), dimension(:,:,:,:), intent(in) :: qn - real(kind=8), dimension(:,:,:,:), intent(out) :: qb + integer , dimension(3) , intent(in) :: pos + real(kind=8), dimension(:,:,:,:), intent(in) :: qn + real(kind=8), dimension(:,:,:,:), intent(inout) :: qb + integer , intent(out) :: status integer :: p - integer :: i, il, iu, is, it, im1, ip1 - integer :: j, jl, ju, js, jt, jm1, jp1 - integer :: k, kt + integer :: i, il, iu, im, ip, is, it + integer :: j, jl, ju, jm, jp, js, jt #if NDIMS == 3 - integer :: kl, ku, ks, km1, kp1 + integer :: k, kl, ku, km, kp, ks, kt +#else /* NDIMS == 3 */ + integer :: k #endif /* NDIMS == 3 */ - real(kind=8) :: dql, dqr real(kind=8) :: dq1, dq2 #if NDIMS == 3 real(kind=8) :: dq3, dq4 @@ -6035,75 +6044,74 @@ module boundaries real(kind=8), dimension(NDIMS) :: dq - character(len=80) :: msg + logical , save :: first = .true. + integer, dimension(2,2), save :: nlims character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong()' -! + !------------------------------------------------------------------------------- ! -#if NDIMS == 2 - il = corners_dp(pos(1),pos(2))%l(1) - jl = corners_dp(pos(1),pos(2))%l(2) - iu = corners_dp(pos(1),pos(2))%u(1) - ju = corners_dp(pos(1),pos(2))%u(2) -#endif /* NDIMS == 2 */ -#if NDIMS == 3 - il = corners_dp(pos(1),pos(2),pos(3))%l(1) - jl = corners_dp(pos(1),pos(2),pos(3))%l(2) - kl = corners_dp(pos(1),pos(2),pos(3))%l(3) - iu = corners_dp(pos(1),pos(2),pos(3))%u(1) - ju = corners_dp(pos(1),pos(2),pos(3))%u(2) - ku = corners_dp(pos(1),pos(2),pos(3))%u(3) -#endif /* NDIMS == 3 */ + if (first) then + nlims(1,1) = ne - nh + 1 + nlims(1,2) = ne + nlims(2,1) = nb + nlims(2,2) = nb + nh - 1 -#if NDIMS == 2 - k = 1 - kt = 1 -#endif /* NDIMS == 2 */ + first = .false. + end if + + il = nlims(pos(1),1) + iu = nlims(pos(1),2) + jl = nlims(pos(2),1) + ju = nlims(pos(2),2) #if NDIMS == 3 + kl = nlims(pos(3),1) + ku = nlims(pos(3),2) + do k = kl, ku - km1 = k - 1 - kp1 = k + 1 - ks = 2 * (k - kl) + 1 - kt = ks + 1 + km = k - 1 + kp = k + 1 + ks = 2 * (k - kl) + 1 + kt = ks + 1 +#else /* NDIMS == 3 */ + k = 1 #endif /* NDIMS == 3 */ do j = jl, ju - jm1 = j - 1 - jp1 = j + 1 - js = 2 * (j - jl) + 1 - jt = js + 1 + jm = j - 1 + jp = j + 1 + js = 2 * (j - jl) + 1 + jt = js + 1 do i = il, iu - im1 = i - 1 - ip1 = i + 1 - is = 2 * (i - il) + 1 - it = is + 1 + im = i - 1 + ip = i + 1 + is = 2 * (i - il) + 1 + it = is + 1 do p = 1, nv - dql = qn(p,i ,j,k) - qn(p,im1,j,k) - dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) - dq(1) = limiter_prol(2.5d-01, dql, dqr) + dq1 = qn(p,i ,j,k) - qn(p,im,j,k) + dq2 = qn(p,ip,j,k) - qn(p,i ,j,k) + dq(1) = limiter_prol(2.5d-01, dq1, dq2) - dql = qn(p,i,j ,k) - qn(p,i,jm1,k) - dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) - dq(2) = limiter_prol(2.5d-01, dql, dqr) + dq1 = qn(p,i,j ,k) - qn(p,i,jm,k) + dq2 = qn(p,i,jp,k) - qn(p,i,j ,k) + dq(2) = limiter_prol(2.5d-01, dq1, dq2) #if NDIMS == 3 - dql = qn(p,i,j,k ) - qn(p,i,j,km1) - dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) - dq(3) = limiter_prol(2.5d-01, dql, dqr) + dq1 = qn(p,i,j,k ) - qn(p,i,j,km) + dq2 = qn(p,i,j,kp) - qn(p,i,j,k ) + dq(3) = limiter_prol(2.5d-01, dq1, dq2) #endif /* NDIMS == 3 */ - if (positive(p) .and. qn(p,i,j,k) < sum(abs(dq(1:NDIMS)))) then + if (positive(p)) then if (qn(p,i,j,k) > 0.0d+00) then - do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) - dq(:) = 0.5d+00 * dq(:) - end do + dq1 = sum(abs(dq(1:NDIMS))) + if (qn(p,i,j,k) <= dq1) & + dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:) else - write(msg,"(a,3i4,a)") & - "Positive variable is not positive at (", i, j, k, " )" - call print_message(loc, msg) - dq(:) = 0.0d+00 + call print_message(loc, "Positive variable is not positive!") + status = 1 + return end if end if @@ -6132,12 +6140,12 @@ module boundaries qb(p,it,jt,kt) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 3 */ - end do ! p + end do - end do ! i - end do ! j + end do + end do #if NDIMS == 3 - end do ! k + end do #endif /* NDIMS == 3 */ !-------------------------------------------------------------------------------