/*@@ @file Staticrad.F @date April 2000 @author Miguel Alcubierre @desc Routine to apply an outgoing radiative boundary condition to the difference between a given function and its initial value. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" subroutine staticrad(bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r, .var_n,var_p,var_0,v0) implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS INTEGER :: i,j,k INTEGER :: nx,ny,nz CCTK_INT, dimension(1:6) :: bbox CCTK_REAL v0,dtv,dtvh CCTK_REAL dt,dx,dy,dz CCTK_REAL rhox,rhoy,rhoz CCTK_REAL, dimension(nx,ny,nz) :: x,y,z,r CCTK_REAL, dimension(nx,ny,nz) :: var_n,var_p,var_0 c Find Courant parameters. dtv = v0*dt dtvh = 0.5D0*dtv rhox = dtv/dx rhoy = dtv/dy rhoz = dtv/dz if (nx.ge.2) then c Lower x. if ((bbox(1).eq.1).and.(.not.CCTK_EQUALS(domain,'octant')) . .and.(.not.CCTK_EQUALS(domain,'quadrant'))) then var_n(2,:,:) = var_n(2,:,:) - var_0(2,:,:) var_p(1,:,:) = var_p(1,:,:) - var_0(1,:,:) var_p(2,:,:) = var_p(2,:,:) - var_0(2,:,:) var_n(1,:,:) = 1.0D0/(rhox - x(1,:,:)/r(1,:,:) . *(1.0D0 + dtvh/r(1,:,:))) . *(var_n(2,:,:)*(rhox + x(2,:,:)/r(2,:,:) . *(1.0D0 + dtvh/r(2,:,:))) . - var_p(1,:,:)*(rhox + x(1,:,:)/r(1,:,:) . *(1.0D0 - dtvh/r(1,:,:))) . + var_p(2,:,:)*(rhox - x(2,:,:)/r(2,:,:) . *(1.0D0 - dtvh/r(2,:,:)))) var_n(1,:,:) = var_n(1,:,:) + var_0(1,:,:) var_n(2,:,:) = var_n(2,:,:) + var_0(2,:,:) var_p(1,:,:) = var_p(1,:,:) + var_0(1,:,:) var_p(2,:,:) = var_p(2,:,:) + var_0(2,:,:) end if c Upper x. if (bbox(2).eq.1) then var_n(nx-1,:,:) = var_n(nx-1,:,:) - var_0(nx-1,:,:) var_p(nx ,:,:) = var_p(nx ,:,:) - var_0(nx ,:,:) var_p(nx-1,:,:) = var_p(nx-1,:,:) - var_0(nx-1,:,:) var_n(nx,:,:) = 1.0D0/(rhox + x(nx ,:,:)/r(nx ,:,:) . *(1.0D0 + dtvh/r(nx ,:,:))) . *(var_n(nx-1,:,:)*(+rhox - x(nx-1,:,:)/r(nx-1,:,:) . *(1.0D0 + dtvh/r(nx-1,:,:))) . + var_p(nx ,:,:)*(-rhox + x(nx ,:,:)/r(nx ,:,:) . *(1.0D0 - dtvh/r(nx ,:,:))) . + var_p(nx-1,:,:)*(+rhox + x(nx-1,:,:)/r(nx-1,:,:) . *(1.0D0 - dtvh/r(nx-1,:,:)))) var_n(nx ,:,:) = var_n(nx ,:,:) + var_0(nx ,:,:) var_n(nx-1,:,:) = var_n(nx-1,:,:) + var_0(nx-1,:,:) var_p(nx ,:,:) = var_p(nx ,:,:) + var_0(nx ,:,:) var_p(nx-1,:,:) = var_p(nx-1,:,:) + var_0(nx-1,:,:) end if end if if (ny.ge.2) then c Lower y. if ((bbox(3).eq.1).and.(.not.CCTK_EQUALS(domain,'octant')) . .and.(.not.CCTK_EQUALS(domain,'quadrant'))) then var_n(:,2,:) = var_n(:,2,:) - var_0(:,2,:) var_p(:,1,:) = var_p(:,1,:) - var_0(:,1,:) var_p(:,2,:) = var_p(:,2,:) - var_0(:,2,:) var_n(:,1,:) = 1.0D0/(rhoy - y(:,1,:)/r(:,1,:) . *(1.0D0 + dtvh/r(:,1,:))) . *(var_n(:,2,:)*(rhoy + y(:,2,:)/r(:,2,:) . *(1.0D0 + dtvh/r(:,2,:))) . - var_p(:,1,:)*(rhoy + y(:,1,:)/r(:,1,:) . *(1.0D0 - dtvh/r(:,1,:))) . + var_p(:,2,:)*(rhoy - y(:,2,:)/r(:,2,:) . *(1.0D0 - dtvh/r(:,2,:)))) var_n(:,1,:) = var_n(:,1,:) + var_0(:,1,:) var_n(:,2,:) = var_n(:,2,:) + var_0(:,2,:) var_p(:,1,:) = var_p(:,1,:) + var_0(:,1,:) var_p(:,2,:) = var_p(:,2,:) + var_0(:,2,:) end if c Upper y. if (bbox(4).eq.1) then var_n(:,ny-1,:) = var_n(:,ny-1,:) - var_0(:,ny-1,:) var_p(:,ny ,:) = var_p(:,ny ,:) - var_0(:,ny ,:) var_p(:,ny-1,:) = var_p(:,ny-1,:) - var_0(:,ny-1,:) var_n(:,ny,:) = 1.0D0/(rhoy + y(:,ny ,:)/r(:,ny ,:) . *(1.0D0 + dtvh/r(:,ny ,:))) . *(var_n(:,ny-1,:)*(+rhoy - y(:,ny-1,:)/r(:,ny-1,:) . *(1.0D0 + dtvh/r(:,ny-1,:))) . + var_p(:,ny ,:)*(-rhoy + y(:,ny ,:)/r(:,ny ,:) . *(1.0D0 - dtvh/r(:,ny ,:))) . + var_p(:,ny-1,:)*(+rhoy + y(:,ny-1,:)/r(:,ny-1,:) . *(1.0D0 - dtvh/r(:,ny-1,:)))) var_n(:,ny ,:) = var_n(:,ny ,:) + var_0(:,ny ,:) var_n(:,ny-1,:) = var_n(:,ny-1,:) + var_0(:,ny-1,:) var_p(:,ny ,:) = var_p(:,ny ,:) + var_0(:,ny ,:) var_p(:,ny-1,:) = var_p(:,ny-1,:) + var_0(:,ny-1,:) end if end if if (nz.ge.2) then c Lower z. if ((bbox(5).eq.1).and.(.not.CCTK_EQUALS(domain,'octant'))) then var_n(:,:,2) = var_n(:,:,2) - var_0(:,:,2) var_p(:,:,1) = var_p(:,:,1) - var_0(:,:,1) var_p(:,:,2) = var_p(:,:,2) - var_0(:,:,2) var_n(:,:,1) = 1.0D0/(rhoz - z(:,:,1)/r(:,:,1) . *(1.0D0 + dtvh/r(:,:,1))) . *(var_n(:,:,2)*(rhoz + z(:,:,2)/r(:,:,2) . *(1.0D0 + dtvh/r(:,:,2))) . - var_p(:,:,1)*(rhoz + z(:,:,1)/r(:,:,1) . *(1.0D0 - dtvh/r(:,:,1))) . + var_p(:,:,2)*(rhoz - z(:,:,2)/r(:,:,2) . *(1.0D0 - dtvh/r(:,:,2)))) var_n(:,:,1) = var_n(:,:,1) + var_0(:,:,1) var_n(:,:,2) = var_n(:,:,2) + var_0(:,:,2) var_p(:,:,1) = var_p(:,:,1) + var_0(:,:,1) var_p(:,:,2) = var_p(:,:,2) + var_0(:,:,2) end if c Upper z. if (bbox(6).eq.1) then var_n(:,:,nz-1) = var_n(:,:,nz-1) - var_0(:,:,nz-1) var_p(:,:,nz ) = var_p(:,:,nz ) - var_0(:,:,nz ) var_p(:,:,nz-1) = var_p(:,:,nz-1) - var_0(:,:,nz-1) var_n(:,:,nz) = 1.0D0/(rhoz + z(:,:,nz )/r(:,:,nz ) . *(1.0D0 + dtvh/r(:,:,nz ))) . *(var_n(:,:,nz-1)*(+rhoz - z(:,:,nz-1)/r(:,:,nz-1) . *(1.0D0 + dtvh/r(:,:,nz-1))) . + var_p(:,:,nz )*(-rhoz + z(:,:,nz )/r(:,:,nz ) . *(1.0D0 - dtvh/r(:,:,nz ))) . + var_p(:,:,nz-1)*(+rhoz + z(:,:,nz-1)/r(:,:,nz-1) . *(1.0D0 - dtvh/r(:,:,nz-1)))) var_n(:,:,nz ) = var_n(:,:,nz ) + var_0(:,:,nz ) var_n(:,:,nz-1) = var_n(:,:,nz-1) + var_0(:,:,nz-1) var_p(:,:,nz ) = var_p(:,:,nz ) + var_0(:,:,nz ) var_p(:,:,nz-1) = var_p(:,:,nz-1) + var_0(:,:,nz-1) end if end if c End. end subroutine Staticrad