From 4634081ca91726ec797cbfc59ab98a135fb0428f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 2 Feb 2022 16:39:51 -0300 Subject: [PATCH] BOUNDARIES: Rewrite boundary_fluxes(). Slighly change the shape of block fluxes. Also add status flag to this subroutine. Signed-off-by: Grzegorz Kowal --- sources/blocks.F90 | 2 +- sources/boundaries.F90 | 1462 +++++++++++++++------------------------- sources/evolution.F90 | 140 +++- 3 files changed, 666 insertions(+), 938 deletions(-) diff --git a/sources/blocks.F90 b/sources/blocks.F90 index 9eb6a96..d11548c 100644 --- a/sources/blocks.F90 +++ b/sources/blocks.F90 @@ -1045,7 +1045,7 @@ module blocks ! allocate(pdata%uu(nvars,nx,ny,nz,nregs), pdata%du(nvars,nx,ny,nz), & pdata%q(nvars,nx,ny,nz), & - pdata%fx(nflux,2,ny,nz), pdata%fy(nflux,nx,2,nz), & + pdata%fx(nflux,ny,nz,2), pdata%fy(nflux,nx,nz,2), & #if NDIMS == 3 pdata%fz(nflux,nx,ny,2), & #endif /* NDIMS == 3 */ diff --git a/sources/boundaries.F90 b/sources/boundaries.F90 index aae2a3f..69e675f 100644 --- a/sources/boundaries.F90 +++ b/sources/boundaries.F90 @@ -297,6 +297,568 @@ module boundaries ! !=============================================================================== ! +! subroutine BOUNDARY_FLUXES: +! -------------------------- +! +! Subroutine scans through the pairs of all neighbors and if one neighbor has +! higher level, takes the boundary flux from it, restricts it to the level of +! the lower level block and corrects its boundary flux. +! +! Arguments: +! +! status - the subroutine call status; +! +!=============================================================================== +! + subroutine boundary_fluxes(status) + + use blocks , only : block_meta, block_data, block_leaf + use blocks , only : list_leaf +#ifdef MPI + use blocks , only : block_info, pointer_info +#endif /* MPI */ + use blocks , only : ndims, nsides + use coordinates, only : minlev, maxlev + use coordinates, only : nh => ncells_half + use coordinates, only : nb, ne, nbm, nbp, nem, nep + use coordinates, only : adxi, adyi +#if NDIMS == 3 + use coordinates, only : adzi +#endif /* NDIMS == 3 */ + use equations , only : nf, ns + use equations , only : idn, isl, isu + use helpers , only : print_message +#ifdef MPI + use mpitools , only : nproc, npairs, pairs + use mpitools , only : exchange_arrays +#endif /* MPI */ + + implicit none + + integer, intent(out) :: status + + type(block_meta), pointer :: pmeta, pneigh + type(block_data), pointer :: pdata + type(block_leaf), pointer :: pleaf +#ifdef MPI + type(block_info), pointer :: pinfo +#endif /* MPI */ + + integer :: n +#if NDIMS == 2 + integer :: m +#endif /* NDIMS == 2 */ + integer :: i, il, iu + integer :: j, jl, ju + integer :: k, kl, ku + integer :: s +#ifdef MPI + integer :: sproc, rproc + integer :: scount, rcount + integer :: l, p +#endif /* MPI */ + +#if NDIMS == 3 + real(kind=8), dimension(nf,nh,nh) :: fh, df + real(kind=8), dimension( nh,nh) :: sl, sr, ps +#else /* NDIMS == 3 */ + real(kind=8), dimension(nf,nh, 1) :: fh, df + real(kind=8), dimension( nh, 1) :: sl, sr, ps +#endif /* NDIMS == 3 */ +#ifdef MPI + real(kind=8), dimension(:,:,:,:), allocatable :: sbuf, rbuf +#endif /* MPI */ + + logical , save :: first = .true. + integer, dimension(2,2), save :: nlims + + character(len=*), parameter :: loc = 'BOUNDARIES::boundary_fluxes()' + +!------------------------------------------------------------------------------- +! + status = 0 + + if (minlev == maxlev) return + + if (first) then + nlims(1,1) = nb + nlims(2,1) = nb + nh + nlims(:,2) = nlims(:,1) + nh - 1 + + first = .false. + end if + +#ifdef MPI + sproc = 0 + rproc = 0 + scount = 0 + rcount = 0 + + call prepare_exchange_array() +#endif /* MPI */ + +#if NDIMS == 2 + k = 1 + kl = 1 + ku = 1 +#endif /* NDIMS == 2 */ + + pleaf => list_leaf + do while(associated(pleaf)) + + pmeta => pleaf%meta + pdata => pmeta%data + + do n = 1, NDIMS +#if NDIMS == 2 + m = 3 - n +#endif /* NDIMS == 2 */ + +#if NDIMS == 3 + do k = 1, nsides +#endif /* NDIMS == 3 */ + do j = 1, nsides + do i = 1, nsides + +#if NDIMS == 2 + pneigh => pmeta%edges(i,j,m)%ptr +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + pneigh => pmeta%faces(i,j,k,n)%ptr +#endif /* NDIMS == 3 */ + + if (associated(pneigh)) then + + if (pneigh%level > pmeta%level) then + +#ifdef MPI + if (pmeta%process == pneigh%process) then + + if (pneigh%process == nproc) then +#endif /* MPI */ + + select case(n) + case(1) + + s = 3 - i + if (i == 1) then + il = nbm + iu = nb + else + il = ne + iu = nep + end if + jl = nlims(j,1) + ju = nlims(j,2) +#if NDIMS == 3 + kl = nlims(k,1) + ku = nlims(k,2) + + fh(:,:,:) = & + 2.5d-01 * ((pneigh%data%fx(:,nb:nem:2,nb:nem:2,s) & + + pneigh%data%fx(:,nbp:ne:2,nbp:ne:2,s)) & + + (pneigh%data%fx(:,nbp:ne:2,nb:nem:2,s) & + + pneigh%data%fx(:,nb:nem:2,nbp:ne:2,s))) +#else /* NDIMS == 3 */ + fh(:,:,:) = & + 5.0d-01 * (pneigh%data%fx(:,nb:nem:2,:,s) & + + pneigh%data%fx(:,nbp:ne:2,:,s)) +#endif /* NDIMS == 3 */ + + df(:,:,:) = (fh(:,:,:) - pdata%fx(:,jl:ju,kl:ku,i)) & + * adxi(pmeta%level) + + pdata%du(:,il,jl:ju,kl:ku) = & + pdata%du(:,il,jl:ju,kl:ku) - df(:,:,:) + pdata%du(:,iu,jl:ju,kl:ku) = & + pdata%du(:,iu,jl:ju,kl:ku) + df(:,:,:) + + if (ns > 0) then + sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 + sr(:,:) = 1.0d+00 - sl(:,:) + do s = isl, isu + ps(:,:) = sl(:,:) * pdata%q(s,il,jl:ju,kl:ku) & + + sr(:,:) * pdata%q(s,iu,jl:ju,kl:ku) + + pdata%du(s,il,jl:ju,kl:ku) = & + pdata%du(s,il,jl:ju,kl:ku) & + - df(idn,:,:) * ps(:,:) + pdata%du(s,iu,jl:ju,kl:ku) = & + pdata%du(s,iu,jl:ju,kl:ku) & + + df(idn,:,:) * ps(:,:) + end do + end if + + case(2) + + s = 3 - j + if (j == 1) then + jl = nbm + ju = nb + else + jl = ne + ju = nep + end if + il = nlims(i,1) + iu = nlims(i,2) +#if NDIMS == 3 + kl = nlims(k,1) + ku = nlims(k,2) + + fh(:,:,:) = & + 2.5d-01 * ((pneigh%data%fy(:,nb:nem:2,nb:nem:2,s) & + + pneigh%data%fy(:,nbp:ne:2,nbp:ne:2,s)) & + + (pneigh%data%fy(:,nbp:ne:2,nb:nem:2,s) & + + pneigh%data%fy(:,nb:nem:2,nbp:ne:2,s))) +#else /* NDIMS == 3 */ + fh(:,:,:) = & + 5.0d-01 * (pneigh%data%fy(:,nb:nem:2,:,s) & + + pneigh%data%fy(:,nbp:ne:2,:,s)) +#endif /* NDIMS == 3 */ + + df(:,:,:) = (fh(:,:,:) - pdata%fy(:,il:iu,kl:ku,j)) & + * adyi(pmeta%level) + + pdata%du(:,il:iu,jl,kl:ku) = & + pdata%du(:,il:iu,jl,kl:ku) - df(:,:,:) + pdata%du(:,il:iu,ju,kl:ku) = & + pdata%du(:,il:iu,ju,kl:ku) + df(:,:,:) + + if (ns > 0) then + sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 + sr(:,:) = 1.0d+00 - sl(:,:) + do s = isl, isu + ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl,kl:ku) & + + sr(:,:) * pdata%q(s,il:iu,ju,kl:ku) + + pdata%du(s,il:iu,jl,kl:ku) = & + pdata%du(s,il:iu,jl,kl:ku) & + - df(idn,:,:) * ps(:,:) + pdata%du(s,il:iu,ju,kl:ku) = & + pdata%du(s,il:iu,ju,kl:ku) & + + df(idn,:,:) * ps(:,:) + end do + end if +#if NDIMS == 3 + + case(3) + + s = 3 - k + if (k == 1) then + kl = nbm + ku = nb + else + kl = ne + ku = nep + end if + il = nlims(i,1) + iu = nlims(i,2) + jl = nlims(j,1) + ju = nlims(j,2) + + fh(:,:,:) = & + 2.5d-01 * ((pneigh%data%fz(:,nb:nem:2,nb:nem:2,s) & + + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,s)) & + + (pneigh%data%fz(:,nbp:ne:2,nb:nem:2,s) & + + pneigh%data%fz(:,nb:nem:2,nbp:ne:2,s))) + + df(:,:,:) = (fh(:,:,:) - pdata%fz(:,il:iu,jl:ju,k)) & + * adzi(pmeta%level) + + pdata%du(:,il:iu,jl:ju,kl) = & + pdata%du(:,il:iu,jl:ju,kl) - df(:,:,:) + pdata%du(:,il:iu,jl:ju,ku) = & + pdata%du(:,il:iu,jl:ju,ku) + df(:,:,:) + + if (ns > 0) then + sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 + sr(:,:) = 1.0d+00 - sl(:,:) + do s = isl, isu + ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,kl) & + + sr(:,:) * pdata%q(s,il:iu,jl:ju,ku) + + pdata%du(s,il:iu,jl:ju,kl) = & + pdata%du(s,il:iu,jl:ju,kl) & + - df(idn,:,:) * ps(:,:) + pdata%du(s,il:iu,jl:ju,ku) = & + pdata%du(s,il:iu,jl:ju,ku) & + + df(idn,:,:) * ps(:,:) + end do + end if +#endif /* NDIMS == 3 */ + end select +#ifdef MPI + end if ! pneigh on the current process + + else ! pneigh%proc /= pmeta%proc + + call append_exchange_block(pmeta, pneigh, n, [ i, j, k ]) + + end if ! pneigh%proc /= pmeta%proc +#endif /* MPI */ + end if ! pmeta level < pneigh level + + end if ! pneigh associated + + end do ! i = 1, nsides + end do ! j = 1, nsides +#if NDIMS == 3 + end do ! k = 1, nsides +#endif /* NDIMS == 3 */ + end do ! n = 1, ndims + + pleaf => pleaf%next + end do ! over leaf blocks + +#ifdef MPI + do p = 1, npairs + + if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then + + if (pairs(p,1) == nproc) then + sproc = pairs(p,1) + rproc = pairs(p,2) + end if + if (pairs(p,2) == nproc) then + sproc = pairs(p,2) + rproc = pairs(p,1) + end if + + scount = bcount(sproc,rproc) + rcount = bcount(rproc,sproc) + + if (scount > 0 .or. rcount > 0) then +#if NDIMS == 3 + allocate(sbuf(nf,nh,nh,scount), rbuf(nf,nh,nh,rcount), stat=status) +#else /* NDIMS == 3 */ + allocate(sbuf(nf,nh, 1,scount), rbuf(nf,nh, 1,rcount), stat=status) +#endif /* NDIMS == 3 */ + if (status /= 0) & + call print_message(loc, "Could not allocate exchange arrays!") + + if (scount > 0) then + + l = 0 + + pinfo => barray(sproc,rproc)%ptr + + do while(associated(pinfo)) + + l = l + 1 + + pdata => pinfo%neigh%data + + select case(pinfo%direction) + case(1) + + s = 3 - pinfo%corner(1) +#if NDIMS == 3 + sbuf(:,:,:,l) = 2.5d-01 * ((pdata%fx(:,nb:nem:2,nb:nem:2,s) & + + pdata%fx(:,nbp:ne:2,nbp:ne:2,s)) & + + (pdata%fx(:,nbp:ne:2,nb:nem:2,s) & + + pdata%fx(:,nb:nem:2,nbp:ne:2,s))) +#else /* NDIMS == 3 */ + sbuf(:,:,:,l) = 5.0d-01 * (pdata%fx(:,nb:nem:2,:,s) & + + pdata%fx(:,nbp:ne:2,:,s)) +#endif /* NDIMS == 3 */ + + case(2) + + s = 3 - pinfo%corner(2) + +#if NDIMS == 3 + sbuf(:,:,:,l) = 2.5d-01 * ((pdata%fy(:,nb:nem:2,nb:nem:2,s) & + + pdata%fy(:,nbp:ne:2,nbp:ne:2,s)) & + + (pdata%fy(:,nbp:ne:2,nb:nem:2,s) & + + pdata%fy(:,nb:nem:2,nbp:ne:2,s))) +#else /* NDIMS == 3 */ + sbuf(:,:,:,l) = 5.0d-01 * (pdata%fy(:,nb:nem:2,:,s) & + + pdata%fy(:,nbp:ne:2,:,s)) +#endif /* NDIMS == 3 */ +#if NDIMS == 3 + + case(3) + + s = 3 - pinfo%corner(3) + + sbuf(:,:,:,l) = 2.5d-01 * ((pdata%fz(:,nb:nem:2,nb:nem:2,s) & + + pdata%fz(:,nbp:ne:2,nbp:ne:2,s)) & + + (pdata%fz(:,nbp:ne:2,nb:nem:2,s) & + + pdata%fz(:,nb:nem:2,nbp:ne:2,s))) +#endif /* NDIMS == 3 */ + + end select + + pinfo => pinfo%prev + end do ! %ptr blocks + + end if + +! exchange_arrays() handles the cases of only send, only receive, or exchange +! of buffers +! + call exchange_arrays(rproc, p, sbuf, rbuf) + + if (rcount > 0) then + + l = 0 + + pinfo => barray(rproc,sproc)%ptr + + do while(associated(pinfo)) + + l = l + 1 + + pmeta => pinfo%meta + pdata => pmeta%data + + n = pinfo%direction + i = pinfo%corner(1) + j = pinfo%corner(2) +#if NDIMS == 3 + k = pinfo%corner(3) +#endif /* NDIMS == 3 */ + + select case(n) + case(1) + + if (i == 1) then + il = nbm + iu = nb + else + il = ne + iu = nep + end if + jl = nlims(j,1) + ju = nlims(j,2) +#if NDIMS == 3 + kl = nlims(k,1) + ku = nlims(k,2) +#endif /* NDIMS == 3 */ + + df(:,:,:) = (rbuf(:,:,:,l) - pdata%fx(:,jl:ju,kl:ku,i)) & + * adxi(pmeta%level) + + pdata%du(:,il,jl:ju,kl:ku) = pdata%du(:,il,jl:ju,kl:ku) & + - df(:,:,:) + pdata%du(:,iu,jl:ju,kl:ku) = pdata%du(:,iu,jl:ju,kl:ku) & + + df(:,:,:) + + if (ns > 0) then + sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 + sr(:,:) = 1.0d+00 - sl(:,:) + do s = isl, isu + ps(:,:) = sl(:,:) * pdata%q(s,il,jl:ju,kl:ku) & + + sr(:,:) * pdata%q(s,iu,jl:ju,kl:ku) + + pdata%du(s,il,jl:ju,kl:ku) = pdata%du(s,il,jl:ju,kl:ku) & + - df(idn,:,:) * ps(:,:) + pdata%du(s,iu,jl:ju,kl:ku) = pdata%du(s,iu,jl:ju,kl:ku) & + + df(idn,:,:) * ps(:,:) + end do + end if + + case(2) + + if (j == 1) then + jl = nbm + ju = nb + else + jl = ne + ju = nep + end if + il = nlims(i,1) + iu = nlims(i,2) +#if NDIMS == 3 + kl = nlims(k,1) + ku = nlims(k,2) +#endif /* NDIMS == 3 */ + + df(:,:,:) = (rbuf(:,:,:,l) - pdata%fy(:,il:iu,kl:ku,j)) & + * adyi(pmeta%level) + + pdata%du(:,il:iu,jl,kl:ku) = pdata%du(:,il:iu,jl,kl:ku) & + - df(:,:,:) + pdata%du(:,il:iu,ju,kl:ku) = pdata%du(:,il:iu,ju,kl:ku) & + + df(:,:,:) + + if (ns > 0) then + sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 + sr(:,:) = 1.0d+00 - sl(:,:) + do s = isl, isu + ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl,kl:ku) & + + sr(:,:) * pdata%q(s,il:iu,ju,kl:ku) + + pdata%du(s,il:iu,jl,kl:ku) = pdata%du(s,il:iu,jl,kl:ku) & + - df(idn,:,:) * ps(:,:) + pdata%du(s,il:iu,ju,kl:ku) = pdata%du(s,il:iu,ju,kl:ku) & + + df(idn,:,:) * ps(:,:) + end do + end if +#if NDIMS == 3 + + case(3) + + if (k == 1) then + kl = nbm + ku = nb + else + kl = ne + ku = nep + end if + il = nlims(i,1) + iu = nlims(i,2) + jl = nlims(j,1) + ju = nlims(j,2) + + df(:,:,:) = (rbuf(:,:,:,l) - pdata%fz(:,il:iu,jl:ju,k)) & + * adzi(pmeta%level) + + pdata%du(:,il:iu,jl:ju,kl) = pdata%du(:,il:iu,jl:ju,kl) & + - df(:,:,:) + pdata%du(:,il:iu,jl:ju,ku) = pdata%du(:,il:iu,jl:ju,ku) & + + df(:,:,:) + + if (ns > 0) then + sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 + sr(:,:) = 1.0d+00 - sl(:,:) + do s = isl, isu + ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,kl) & + + sr(:,:) * pdata%q(s,il:iu,jl:ju,ku) + + pdata%du(s,il:iu,jl:ju,kl) = pdata%du(s,il:iu,jl:ju,kl) & + - df(idn,:,:) * ps(:,:) + pdata%du(s,il:iu,jl:ju,ku) = pdata%du(s,il:iu,jl:ju,ku) & + + df(idn,:,:) * ps(:,:) + end do + end if +#endif /* NDIMS == 3 */ + + end select + + pinfo => pinfo%prev + end do ! %ptr blocks + end if + + deallocate(sbuf, rbuf, stat=status) + if (status /= 0) & + call print_message(loc, "Could not deallocate exchange arrays!") + + end if + + end if ! pairs(p,1) == nproc || pairs(p,2) == nproc + + end do ! p = 1, npairs + + call release_exchange_array() +#endif /* MPI */ + +!------------------------------------------------------------------------------- +! + end subroutine boundary_fluxes +! +!=============================================================================== +! ! subroutine BOUNDARY_VARIABLES: ! ----------------------------- ! @@ -401,906 +963,6 @@ module boundaries end subroutine boundary_variables ! !=============================================================================== -! -! subroutine BOUNDARY_FLUXES: -! -------------------------- -! -! Subroutine updates the numerical fluxes from neighbors which lay on -! higher level. At higher levels the numerical fluxes are calculated with -! smaller error, since the resolution is higher, therefore we take those -! fluxes and restrict them down to the level of the updated block. -! -! -!=============================================================================== -! - subroutine boundary_fluxes() - -! import external procedures and variables -! - use blocks , only : block_meta, block_data, block_leaf - use blocks , only : list_leaf -#ifdef MPI - use blocks , only : block_info, pointer_info -#endif /* MPI */ - use blocks , only : ndims, nsides - use coordinates, only : minlev, maxlev - use coordinates, only : nh => ncells_half - use coordinates, only : nb, ne, nbm, nbp, nep - use coordinates, only : adxi, adyi -#if NDIMS == 3 - use coordinates, only : adzi -#endif /* NDIMS == 3 */ - use equations , only : nf, ns - use equations , only : idn, isl, isu -#ifdef MPI - use mpitools , only : nproc, npairs, pairs - use mpitools , only : exchange_arrays -#endif /* MPI */ - -! local variables are not implicit by default -! - implicit none - -! local pointers -! - type(block_meta), pointer :: pmeta, pneigh - type(block_data), pointer :: pdata - type(block_leaf), pointer :: pleaf -#ifdef MPI - type(block_info), pointer :: pinfo -#endif /* MPI */ - -! local variables -! - integer :: n -#if NDIMS == 2 - integer :: m -#endif /* NDIMS == 2 */ - integer :: i, il, iu - integer :: j, jl, ju - integer :: k, kl, ku - integer :: s -#ifdef MPI - integer :: sproc, rproc - integer :: scount, rcount - integer :: l, p - -! local arrays -! - real(kind=8), dimension(:,:,:,:), allocatable :: sbuf, rbuf -#endif /* MPI */ -#if NDIMS == 3 - real(kind=8), dimension(nf,nh,nh) :: fh, df - real(kind=8), dimension(nh,nh) :: sl, sr, ps -#else /* NDIMS == 3 */ - real(kind=8), dimension(nf,nh, 1) :: fh, df - real(kind=8), dimension(nh, 1) :: sl, sr, ps -#endif /* NDIMS == 3 */ -! -!------------------------------------------------------------------------------- -! -! quit if all blocks are at the same level -! - if (minlev == maxlev) return - -#ifdef MPI - sproc = 0 - rproc = 0 - scount = 0 - rcount = 0 - -! prepare the block exchange structures -! - call prepare_exchange_array() -#endif /* MPI */ - -#if NDIMS == 2 - k = 1 - kl = 1 - ku = 1 -#endif /* NDIMS == 2 */ - -! update the fluxes between blocks on the same process -! -! associate pleaf with the first block on the leaf list -! - pleaf => list_leaf - -! scan all leaf meta blocks in the list -! - do while(associated(pleaf)) - -! get the associated meta and data block -! - pmeta => pleaf%meta - pdata => pmeta%data - -! iterate over all dimensions -! - do n = 1, ndims -#if NDIMS == 2 - m = 3 - n -#endif /* NDIMS == 2 */ - -! iterate over all corners -! -#if NDIMS == 3 - do k = 1, nsides -#endif /* NDIMS == 3 */ - do j = 1, nsides - do i = 1, nsides - -! associate pneigh with the neighbor -! -#if NDIMS == 2 - pneigh => pmeta%edges(i,j,m)%ptr -#endif /* NDIMS == 2 */ -#if NDIMS == 3 - pneigh => pmeta%faces(i,j,k,n)%ptr -#endif /* NDIMS == 3 */ - -! process only if the neighbor is associated -! - if (associated(pneigh)) then - -! check if the neighbor lays at higher level -! - if (pneigh%level > pmeta%level) then - -#ifdef MPI -! check if the block and its neighbor belong to the same process -! - if (pmeta%process == pneigh%process) then - -! check if the neighbor belongs to the current process -! - if (pneigh%process == nproc) then -#endif /* MPI */ - -! update the flux depending on the direction -! - select case(n) - case(1) - -! prepare the boundary layer indices for X-direction flux -! - if (j == 1) then - jl = nb - ju = nb + nh - 1 - else - jl = ne - nh + 1 - ju = ne - end if -#if NDIMS == 3 - if (k == 1) then - kl = nb - ku = nb + nh - 1 - else - kl = ne - nh + 1 - ku = ne - end if -#endif /* NDIMS == 3 */ - -! update the flux at the X-face of the block -! - if (i == 1) then -#if NDIMS == 3 - fh(:,:,:) = & - 2.5d-01 * ((pneigh%data%fx(:,2,nb :ne:2,nb :ne:2) & - + pneigh%data%fx(:,2,nbp:ne:2,nbp:ne:2)) & - + (pneigh%data%fx(:,2,nbp:ne:2,nb :ne:2) & - + pneigh%data%fx(:,2,nb :ne:2,nbp:ne:2))) -#else /* NDIMS == 3 */ - fh(:,:,:) = & - 5.0d-01 * (pneigh%data%fx(:,2,nb :ne:2,:) & - + pneigh%data%fx(:,2,nbp:ne:2,:)) -#endif /* NDIMS == 3 */ - df(:,:,:) = (fh(:,:,:) - pdata%fx(:,1,jl:ju,kl:ku)) & - * adxi(pmeta%level) - - pdata%du(:,nbm,jl:ju,kl:ku) = & - pdata%du(:,nbm,jl:ju,kl:ku) - df(:,:,:) - pdata%du(:,nb ,jl:ju,kl:ku) = & - pdata%du(:,nb ,jl:ju,kl:ku) + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,nbm,jl:ju,kl:ku) & - + sr(:,:) * pdata%q(s,nb ,jl:ju,kl:ku) - - pdata%du(s,nbm,jl:ju,kl:ku) = & - pdata%du(s,nbm,jl:ju,kl:ku) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,nb ,jl:ju,kl:ku) = & - pdata%du(s,nb ,jl:ju,kl:ku) & - + df(idn,:,:) * ps(:,:) - end do - end if - else -#if NDIMS == 3 - fh(:,:,:) = & - 2.5d-01 * ((pneigh%data%fx(:,1,nb :ne:2,nb :ne:2) & - + pneigh%data%fx(:,1,nbp:ne:2,nbp:ne:2)) & - + (pneigh%data%fx(:,1,nbp:ne:2,nb :ne:2) & - + pneigh%data%fx(:,1,nb :ne:2,nbp:ne:2))) -#else /* NDIMS == 3 */ - fh(:,:,:) = & - 5.0d-01 * (pneigh%data%fx(:,1,nb :ne:2,:) & - + pneigh%data%fx(:,1,nbp:ne:2,:)) -#endif /* NDIMS == 3 */ - df(:,:,:) = (fh(:,:,:) - pdata%fx(:,2,jl:ju,kl:ku)) & - * adxi(pmeta%level) - - pdata%du(:,ne ,jl:ju,kl:ku) = & - pdata%du(:,ne ,jl:ju,kl:ku) - df(:,:,:) - pdata%du(:,nep,jl:ju,kl:ku) = & - pdata%du(:,nep,jl:ju,kl:ku) + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,ne ,jl:ju,kl:ku) & - + sr(:,:) * pdata%q(s,nep,jl:ju,kl:ku) - - pdata%du(s,ne ,jl:ju,kl:ku) = & - pdata%du(s,ne ,jl:ju,kl:ku) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,nep,jl:ju,kl:ku) = & - pdata%du(s,nep,jl:ju,kl:ku) & - + df(idn,:,:) * ps(:,:) - end do - end if - end if - - case(2) - -! prepare the boundary layer indices for Y-direction flux -! - if (i == 1) then - il = nb - iu = nb + nh - 1 - else - il = ne - nh + 1 - iu = ne - end if -#if NDIMS == 3 - if (k == 1) then - kl = nb - ku = nb + nh - 1 - else - kl = ne - nh + 1 - ku = ne - end if -#endif /* NDIMS == 3 */ - -! update the flux at the Y-face of the block -! - if (j == 1) then -#if NDIMS == 3 - fh(:,:,:) = & - 2.5d-01 * ((pneigh%data%fy(:,nb :ne:2,2,nb :ne:2) & - + pneigh%data%fy(:,nbp:ne:2,2,nbp:ne:2)) & - + (pneigh%data%fy(:,nbp:ne:2,2,nb :ne:2) & - + pneigh%data%fy(:,nb :ne:2,2,nbp:ne:2))) -#else /* NDIMS == 3 */ - fh(:,:,:) = & - 5.0d-01 * (pneigh%data%fy(:,nb :ne:2,2,:) & - + pneigh%data%fy(:,nbp:ne:2,2,:)) -#endif /* NDIMS == 3 */ - df(:,:,:) = (fh(:,:,:) - pdata%fy(:,il:iu,1,kl:ku)) & - * adyi(pmeta%level) - - pdata%du(:,il:iu,nbm,kl:ku) = & - pdata%du(:,il:iu,nbm,kl:ku) - df(:,:,:) - pdata%du(:,il:iu,nb ,kl:ku) = & - pdata%du(:,il:iu,nb ,kl:ku) + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,il:iu,nbm,kl:ku) & - + sr(:,:) * pdata%q(s,il:iu,nb ,kl:ku) - - pdata%du(s,il:iu,nbm,kl:ku) = & - pdata%du(s,il:iu,nbm,kl:ku) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,il:iu,nb ,kl:ku) = & - pdata%du(s,il:iu,nb ,kl:ku) & - + df(idn,:,:) * ps(:,:) - end do - end if - else -#if NDIMS == 3 - fh(:,:,:) = & - 2.5d-01 * ((pneigh%data%fy(:,nb :ne:2,1,nb :ne:2) & - + pneigh%data%fy(:,nbp:ne:2,1,nbp:ne:2)) & - + (pneigh%data%fy(:,nbp:ne:2,1,nb :ne:2) & - + pneigh%data%fy(:,nb :ne:2,1,nbp:ne:2))) -#else /* NDIMS == 3 */ - fh(:,:,:) = & - 5.0d-01 * (pneigh%data%fy(:,nb :ne:2,1,:) & - + pneigh%data%fy(:,nbp:ne:2,1,:)) -#endif /* NDIMS == 3 */ - df(:,:,:) = (fh(:,:,:) - pdata%fy(:,il:iu,2,kl:ku)) & - * adyi(pmeta%level) - - pdata%du(:,il:iu,ne ,kl:ku) = & - pdata%du(:,il:iu,ne ,kl:ku) - df(:,:,:) - pdata%du(:,il:iu,nep,kl:ku) = & - pdata%du(:,il:iu,nep,kl:ku) + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,il:iu,ne ,kl:ku) & - + sr(:,:) * pdata%q(s,il:iu,nep,kl:ku) - - pdata%du(s,il:iu,ne ,kl:ku) = & - pdata%du(s,il:iu,ne ,kl:ku) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,il:iu,nep,kl:ku) = & - pdata%du(s,il:iu,nep,kl:ku) & - + df(idn,:,:) * ps(:,:) - end do - end if - end if - -#if NDIMS == 3 - case(3) - -! prepare the boundary layer indices for Z-direction flux -! - if (i == 1) then - il = nb - iu = nb + nh - 1 - else - il = ne - nh + 1 - iu = ne - end if - if (j == 1) then - jl = nb - ju = nb + nh - 1 - else - jl = ne - nh + 1 - ju = ne - end if - -! update the flux at the Z-face of the block -! - if (k == 1) then - fh(:,:,:) = & - 2.5d-01 * ((pneigh%data%fz(:,nb :ne:2,nb :ne:2,2) & - + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,2)) & - + (pneigh%data%fz(:,nbp:ne:2,nb :ne:2,2) & - + pneigh%data%fz(:,nb :ne:2,nbp:ne:2,2))) - df(:,:,:) = (fh(:,:,:) - pdata%fz(:,il:iu,jl:ju,1)) & - * adzi(pmeta%level) - - pdata%du(:,il:iu,jl:ju,nbm) = & - pdata%du(:,il:iu,jl:ju,nbm) - df(:,:,:) - pdata%du(:,il:iu,jl:ju,nb ) = & - pdata%du(:,il:iu,jl:ju,nb ) + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,nbm) & - + sr(:,:) * pdata%q(s,il:iu,jl:ju,nb ) - - pdata%du(s,il:iu,jl:ju,nbm) = & - pdata%du(s,il:iu,jl:ju,nbm) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,il:iu,jl:ju,nb ) = & - pdata%du(s,il:iu,jl:ju,nb ) & - + df(idn,:,:) * ps(:,:) - end do - end if - else - fh(:,:,:) = & - 2.5d-01 * ((pneigh%data%fz(:,nb :ne:2,nb :ne:2,1) & - + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,1)) & - + (pneigh%data%fz(:,nbp:ne:2,nb :ne:2,1) & - + pneigh%data%fz(:,nb :ne:2,nbp:ne:2,1))) - df(:,:,:) = (fh(:,:,:) - pdata%fz(:,il:iu,jl:ju,2)) & - * adzi(pmeta%level) - - pdata%du(:,il:iu,jl:ju,ne ) = & - pdata%du(:,il:iu,jl:ju,ne ) - df(:,:,:) - pdata%du(:,il:iu,jl:ju,nep) = & - pdata%du(:,il:iu,jl:ju,nep) + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,ne ) & - + sr(:,:) * pdata%q(s,il:iu,jl:ju,nep) - - pdata%du(s,il:iu,jl:ju,ne ) = & - pdata%du(s,il:iu,jl:ju,ne ) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,il:iu,jl:ju,nep) = & - pdata%du(s,il:iu,jl:ju,nep) & - + df(idn,:,:) * ps(:,:) - end do - end if - end if -#endif /* NDIMS == 3 */ - - end select - -#ifdef MPI - end if ! pneigh on the current process - - else ! pneigh%proc /= pmeta%proc - -! append the block to the exchange list -! - call append_exchange_block(pmeta, pneigh, n, (/ i, j, k /)) - - end if ! pneigh%proc /= pmeta%proc -#endif /* MPI */ - end if ! pmeta level < pneigh level - - end if ! pneigh associated - - end do ! i = 1, nsides - end do ! j = 1, nsides -#if NDIMS == 3 - end do ! k = 1, nsides -#endif /* NDIMS == 3 */ - end do ! n = 1, ndims - -! associate pleaf with the next leaf on the list -! - pleaf => pleaf%next - - end do ! over leaf blocks - -#ifdef MPI -! update flux boundaries between neighbors laying on different processes -! -! iterate over all process pairs -! - do p = 1, npairs - -! process only pairs related to this process -! - if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then - -! get sending and receiving process identifiers (depending on pair member) -! - if (pairs(p,1) == nproc) then - sproc = pairs(p,1) - rproc = pairs(p,2) - end if - if (pairs(p,2) == nproc) then - sproc = pairs(p,2) - rproc = pairs(p,1) - end if - -! get the number of blocks to exchange -! - scount = bcount(sproc,rproc) - rcount = bcount(rproc,sproc) - -! process only pairs which have anything to exchange -! - if ((scount + rcount) > 0) then - -! allocate buffers for variable exchange -! -#if NDIMS == 3 - allocate(sbuf(nf,nh,nh,scount)) - allocate(rbuf(nf,nh,nh,rcount)) -#else /* NDIMS == 3 */ - allocate(sbuf(nf,nh, 1,scount)) - allocate(rbuf(nf,nh, 1,rcount)) -#endif /* NDIMS == 3 */ - -!! PREPARE BLOCKS FOR SENDING -!! -! reset the block counter -! - l = 0 - -! associate pinfo with the first block in the exchange list -! - pinfo => barray(sproc,rproc)%ptr - -! scan all blocks on the list -! - do while(associated(pinfo)) - -! increase the block count -! - l = l + 1 - -! associate pneigh pointer -! - pneigh => pinfo%neigh - -! get neighbor direction and corner coordinates -! - n = pinfo%direction - i = pinfo%corner(1) - j = pinfo%corner(2) -#if NDIMS == 3 - k = pinfo%corner(3) -#endif /* NDIMS == 3 */ - -! update directional flux from the neighbor -! - select case(n) - case(1) - -! update the flux edge from the neighbor at higher level -! - if (i == 1) then -#if NDIMS == 3 - sbuf(:,:,:,l) = & - 2.5d-01 * ((pneigh%data%fx(:,2,nb :ne:2,nb :ne:2) & - + pneigh%data%fx(:,2,nbp:ne:2,nbp:ne:2)) & - + (pneigh%data%fx(:,2,nbp:ne:2,nb :ne:2) & - + pneigh%data%fx(:,2,nb :ne:2,nbp:ne:2))) -#else /* NDIMS == 3 */ - sbuf(:,:,:,l) = & - 5.0d-01 * (pneigh%data%fx(:,2,nb :ne:2,:) & - + pneigh%data%fx(:,2,nbp:ne:2,:)) -#endif /* NDIMS == 3 */ - else -#if NDIMS == 3 - sbuf(:,:,:,l) = & - 2.5d-01 * ((pneigh%data%fx(:,1,nb :ne:2,nb :ne:2) & - + pneigh%data%fx(:,1,nbp:ne:2,nbp:ne:2)) & - + (pneigh%data%fx(:,1,nbp:ne:2,nb :ne:2) & - + pneigh%data%fx(:,1,nb :ne:2,nbp:ne:2))) -#else /* NDIMS == 3 */ - sbuf(:,:,:,l) = & - 5.0d-01 * (pneigh%data%fx(:,1,nb :ne:2,:) & - + pneigh%data%fx(:,1,nbp:ne:2,:)) -#endif /* NDIMS == 3 */ - end if - - case(2) - -! update the flux edge from the neighbor at higher level -! - if (j == 1) then -#if NDIMS == 3 - sbuf(:,:,:,l) = & - 2.5d-01 * ((pneigh%data%fy(:,nb :ne:2,2,nb :ne:2) & - + pneigh%data%fy(:,nbp:ne:2,2,nbp:ne:2)) & - + (pneigh%data%fy(:,nbp:ne:2,2,nb :ne:2) & - + pneigh%data%fy(:,nb :ne:2,2,nbp:ne:2))) -#else /* NDIMS == 3 */ - sbuf(:,:,:,l) = & - 5.0d-01 * (pneigh%data%fy(:,nb :ne:2,2,:) & - + pneigh%data%fy(:,nbp:ne:2,2,:)) -#endif /* NDIMS == 3 */ - else -#if NDIMS == 3 - sbuf(:,:,:,l) = & - 2.5d-01 * ((pneigh%data%fy(:,nb :ne:2,1,nb :ne:2) & - + pneigh%data%fy(:,nbp:ne:2,1,nbp:ne:2)) & - + (pneigh%data%fy(:,nbp:ne:2,1,nb :ne:2) & - + pneigh%data%fy(:,nb :ne:2,1,nbp:ne:2))) -#else /* NDIMS == 3 */ - sbuf(:,:,:,l) = & - 5.0d-01 * (pneigh%data%fy(:,nb :ne:2,1,:) & - + pneigh%data%fy(:,nbp:ne:2,1,:)) -#endif /* NDIMS == 3 */ - end if - -#if NDIMS == 3 - case(3) - -! update the flux edge from the neighbor at higher level -! - if (k == 1) then - sbuf(:,:,:,l) = & - 2.5d-01 * ((pneigh%data%fz(:,nb :ne:2,nb :ne:2,2) & - + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,2)) & - + (pneigh%data%fz(:,nbp:ne:2,nb :ne:2,2) & - + pneigh%data%fz(:,nb :ne:2,nbp:ne:2,2))) - else - sbuf(:,:,:,l) = & - 2.5d-01 * ((pneigh%data%fz(:,nb :ne:2,nb :ne:2,1) & - + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,1)) & - + (pneigh%data%fz(:,nbp:ne:2,nb :ne:2,1) & - + pneigh%data%fz(:,nb :ne:2,nbp:ne:2,1))) - end if -#endif /* NDIMS == 3 */ - - end select - -! associate pinfo with the next block -! - pinfo => pinfo%prev - - end do ! %ptr blocks - -!! SEND PREPARED BLOCKS AND RECEIVE NEW ONES -!! -! exchange data -! - call exchange_arrays(rproc, p, sbuf, rbuf) - -!! PROCESS RECEIVED BLOCKS -!! -! reset the block counter -! - l = 0 - -! associate pinfo with the first block in the exchange list -! - pinfo => barray(rproc,sproc)%ptr - -! scan all blocks on the list -! - do while(associated(pinfo)) - -! increase the block count -! - l = l + 1 - -! associate meta and data block pointers -! - pmeta => pinfo%meta - pdata => pmeta%data - -! get neighbor direction and corner indices -! - n = pinfo%direction - i = pinfo%corner(1) - j = pinfo%corner(2) -#if NDIMS == 3 - k = pinfo%corner(3) -#endif /* NDIMS == 3 */ - -! update directional flux from the neighbor -! - select case(n) - case(1) - -! prepare the boundary layer indices depending on the corner position -! - if (j == 1) then - jl = nb - ju = nb + nh - 1 - else - jl = ne - nh + 1 - ju = ne - end if -#if NDIMS == 3 - if (k == 1) then - kl = nb - ku = nb + nh - 1 - else - kl = ne - nh + 1 - ku = ne - end if -#endif /* NDIMS == 3 */ - -! update the flux edge from the neighbor at higher level -! - if (i == 1) then - df(:,:,:) = (rbuf(:,:,:,l) - pdata%fx(:,1,jl:ju,kl:ku)) & - * adxi(pmeta%level) - - pdata%du(:,nbm,jl:ju,kl:ku) = & - pdata%du(:,nbm,jl:ju,kl:ku) - df(:,:,:) - pdata%du(:,nb ,jl:ju,kl:ku) = & - pdata%du(:,nb ,jl:ju,kl:ku) + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,nbm,jl:ju,kl:ku) & - + sr(:,:) * pdata%q(s,nb ,jl:ju,kl:ku) - - pdata%du(s,nbm,jl:ju,kl:ku) = pdata%du(s,nbm,jl:ju,kl:ku) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,nb ,jl:ju,kl:ku) = pdata%du(s,nb ,jl:ju,kl:ku) & - + df(idn,:,:) * ps(:,:) - end do - end if - else - df(:,:,:) = (rbuf(:,:,:,l) - pdata%fx(:,2,jl:ju,kl:ku)) & - * adxi(pmeta%level) - - pdata%du(:,ne ,jl:ju,kl:ku) = & - pdata%du(:,ne ,jl:ju,kl:ku) - df(:,:,:) - pdata%du(:,nep,jl:ju,kl:ku) = & - pdata%du(:,nep,jl:ju,kl:ku) + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,ne ,jl:ju,kl:ku) & - + sr(:,:) * pdata%q(s,nep,jl:ju,kl:ku) - - pdata%du(s,ne ,jl:ju,kl:ku) = pdata%du(s,ne ,jl:ju,kl:ku) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,nep,jl:ju,kl:ku) = pdata%du(s,nep,jl:ju,kl:ku) & - + df(idn,:,:) * ps(:,:) - end do - end if - end if - - case(2) - -! prepare the boundary layer indices depending on the corner position -! - if (i == 1) then - il = nb - iu = nb + nh - 1 - else - il = ne - nh + 1 - iu = ne - end if -#if NDIMS == 3 - if (k == 1) then - kl = nb - ku = nb + nh - 1 - else - kl = ne - nh + 1 - ku = ne - end if -#endif /* NDIMS == 3 */ - -! update the flux edge from the neighbor at higher level -! - if (j == 1) then - df(:,:,:) = (rbuf(:,:,:,l) - pdata%fy(:,il:iu,1,kl:ku)) & - * adyi(pmeta%level) - - pdata%du(:,il:iu,nbm,kl:ku) = pdata%du(:,il:iu,nbm,kl:ku) & - - df(:,:,:) - pdata%du(:,il:iu,nb ,kl:ku) = pdata%du(:,il:iu,nb ,kl:ku) & - + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,il:iu,nbm,kl:ku) & - + sr(:,:) * pdata%q(s,il:iu,nb ,kl:ku) - - pdata%du(s,il:iu,nbm,kl:ku) = pdata%du(s,il:iu,nbm,kl:ku) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,il:iu,nb ,kl:ku) = pdata%du(s,il:iu,nb ,kl:ku) & - + df(idn,:,:) * ps(:,:) - end do - end if - else - df(:,:,:) = (rbuf(:,:,:,l) - pdata%fy(:,il:iu,2,kl:ku)) & - * adyi(pmeta%level) - - pdata%du(:,il:iu,ne ,kl:ku) = pdata%du(:,il:iu,ne ,kl:ku) & - - df(:,:,:) - pdata%du(:,il:iu,nep,kl:ku) = pdata%du(:,il:iu,nep,kl:ku) & - + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,il:iu,ne ,kl:ku) & - + sr(:,:) * pdata%q(s,il:iu,nep,kl:ku) - - pdata%du(s,il:iu,ne ,kl:ku) = pdata%du(s,il:iu,ne ,kl:ku) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,il:iu,nep,kl:ku) = pdata%du(s,il:iu,nep,kl:ku) & - + df(idn,:,:) * ps(:,:) - end do - end if - end if - -#if NDIMS == 3 - case(3) - -! prepare the boundary layer indices depending on the corner position -! - if (i == 1) then - il = nb - iu = nb + nh - 1 - else - il = ne - nh + 1 - iu = ne - end if - if (j == 1) then - jl = nb - ju = nb + nh - 1 - else - jl = ne - nh + 1 - ju = ne - end if - -! update the flux edge from the neighbor at higher level -! - if (k == 1) then - df(:,:,:) = (rbuf(:,:,:,l) - pdata%fz(:,il:iu,jl:ju,1)) & - * adzi(pmeta%level) - - pdata%du(:,il:iu,jl:ju,nbm) = pdata%du(:,il:iu,jl:ju,nbm) & - - df(:,:,:) - pdata%du(:,il:iu,jl:ju,nb ) = pdata%du(:,il:iu,jl:ju,nb ) & - + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,nbm) & - + sr(:,:) * pdata%q(s,il:iu,jl:ju,nb ) - - pdata%du(s,il:iu,jl:ju,nbm) = pdata%du(s,il:iu,jl:ju,nbm) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,il:iu,jl:ju,nb ) = pdata%du(s,il:iu,jl:ju,nb ) & - + df(idn,:,:) * ps(:,:) - end do - end if - else - df(:,:,:) = (rbuf(:,:,:,l) - pdata%fz(:,il:iu,jl:ju,2)) & - * adzi(pmeta%level) - pdata%du(:,il:iu,jl:ju,ne ) = pdata%du(:,il:iu,jl:ju,ne ) & - - df(:,:,:) - pdata%du(:,il:iu,jl:ju,nep) = pdata%du(:,il:iu,jl:ju,nep) & - + df(:,:,:) - - if (ns > 0) then - sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 - sr(:,:) = 1.0d+00 - sl(:,:) - do s = isl, isu - ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,ne ) & - + sr(:,:) * pdata%q(s,il:iu,jl:ju,nep) - - pdata%du(s,il:iu,jl:ju,ne ) = pdata%du(s,il:iu,jl:ju,ne ) & - - df(idn,:,:) * ps(:,:) - pdata%du(s,il:iu,jl:ju,nep) = pdata%du(s,il:iu,jl:ju,nep) & - + df(idn,:,:) * ps(:,:) - end do - end if - end if -#endif /* NDIMS == 3 */ - - end select - -! associate pinfo with the next block -! - pinfo => pinfo%prev - - end do ! %ptr blocks - -! deallocate data buffer -! - deallocate(sbuf, rbuf) - - end if ! (scount + rcount) > 0 - - end if ! pairs(p,1) == nproc || pairs(p,2) == nproc - - end do ! p = 1, npairs - -! release the memory used by the array of exchange block lists -! - call release_exchange_array() -#endif /* MPI */ - -!------------------------------------------------------------------------------- -! - end subroutine boundary_fluxes -! -!=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 9d171de..8d50122 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1530,7 +1530,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) return !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1625,7 +1627,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1652,7 +1656,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1752,7 +1758,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1782,7 +1790,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1811,7 +1821,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1915,7 +1927,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1944,7 +1958,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1971,7 +1987,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2000,7 +2018,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2109,7 +2129,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2138,7 +2160,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2166,7 +2190,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2195,7 +2221,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2226,7 +2254,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2347,7 +2377,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2377,7 +2409,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2407,7 +2441,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2438,7 +2474,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2471,7 +2509,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2594,7 +2634,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2640,7 +2682,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2669,7 +2713,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2806,7 +2852,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -2839,7 +2887,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -3038,7 +3088,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -3083,7 +3135,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -3308,7 +3362,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -3349,7 +3405,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -3400,7 +3458,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -3593,7 +3653,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 end if !$omp parallel do default(shared) private(pdata) @@ -3627,7 +3689,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -3682,7 +3746,9 @@ module evolution end do !$omp end parallel do - call boundary_fluxes() + call boundary_fluxes(status) + + if (status /= 0) go to 100 !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -3876,10 +3942,10 @@ module evolution ! store the block interface fluxes ! - pdata%fx(:,1,:,:) = f(:,nbl,:,:,1) - pdata%fx(:,2,:,:) = f(:,ne ,:,:,1) - pdata%fy(:,:,1,:) = f(:,:,nbl,:,2) - pdata%fy(:,:,2,:) = f(:,:,ne ,:,2) + pdata%fx(:,:,:,1) = f(:,nbl,:,:,1) + pdata%fx(:,:,:,2) = f(:,ne ,:,:,1) + pdata%fy(:,:,:,1) = f(:,:,nbl,:,2) + pdata%fy(:,:,:,2) = f(:,:,ne ,:,2) #if NDIMS == 3 pdata%fz(:,:,:,1) = f(:,:,:,nbl,3) pdata%fz(:,:,:,2) = f(:,:,:,ne ,3)