! ! Correct fluxes at boundaries for subsonic flows to avoid ! ill posedness and allow large temporal integrations. ! Creator: Carlos Pantano ! Copyright: California Institute of Technology (2004-2006) ! Notes: This module is part of the CLES patch solver. MODULE Generic_AcousticBC SAVE INTEGER CLESLOG_CBC INTEGER CLESLOG_CBC_EIGENSYSTEM INTEGER CLESLOG_CBC_REFPRESSURE PARAMETER (CLESLOG_CBC=1) PARAMETER (CLESLOG_CBC_EIGENSYSTEM=1) PARAMETER (CLESLOG_CBC_REFPRESSURE=2) ! Determine if acoustic BC should be use ! Use characteristic acoustic BC and incomming as suggested by ! Thompson; Rudy & Strikwerda; Poinsot & Lele; Stanley & Sarkar; ! Giles; Rowley & Colonius. (many theoretical contributions) but ! using the SAT procesure of Carpenter, Gotlieb and Abarbanel. INTEGER :: cbc_active = 0 INTEGER :: cbc_mode = 2 ! default is 2d mode DOUBLE PRECISION, ALLOCATABLE :: cbc_force_matrix(:,:) DOUBLE PRECISION :: cbc_sigma = 0.25d0 DOUBLE PRECISION :: cbc_Mach = 0.0d0 DOUBLE PRECISION :: cbc_factor = 0.125d0 ! ---- Adds acoustic forcing to the state flux to mimmic the ! ---- physical effects that a real infinite domain would ! ---- have on a subsonic outflow. INTERFACE AcousticBC MODULE PROCEDURE OneDAcousticBC MODULE PROCEDURE TwoDAcousticBC MODULE PROCEDURE ThreeDAcousticBC END INTERFACE INTERFACE AcousticViscousBC MODULE PROCEDURE OneDAcousticViscousBC MODULE PROCEDURE TwoDAcousticViscousBC MODULE PROCEDURE ThreeDAcousticViscousBC END INTERFACE CONTAINS SUBROUTINE SingleCSAT(ux, vx, uxref, source, sat, flux, p, dl, tau, & cbc_type, ipos, jpos, kpos) ! ---- Shared Variables USE mesh USE array_bounds USE method_parms use cles_interfaces USE Generic_Source USE Generic_EigenSystem USE Generic_OuterBoundary USE Interp_coeffs IMPLICIT NONE include 'cles.i' DOUBLE PRECISION, INTENT(IN) :: ux(ncomps), dl, p DOUBLE PRECISION, INTENT(IN) :: vx(nvars), tau(nvars) DOUBLE PRECISION, INTENT(IN) :: uxref(ncomps) DOUBLE PRECISION, INTENT(INOUT) :: flux(nvars,3) DOUBLE PRECISION, INTENT(INOUT) :: source(nvars) DOUBLE PRECISION, INTENT(OUT) :: sat(nvars) INTEGER, INTENT(IN) :: ipos, jpos, kpos, cbc_type DOUBLE PRECISION :: a, c, fx, fy, fz DOUBLE PRECISION :: lamda(nvars), Rdir(nvars,nvars), Rinv(nvars,nvars) DOUBLE PRECISION :: vec(nvars), dU(nvars), DeltaG(nvars) INTEGER :: i, k CHARACTER(len=56) :: slog call cleslog_log_enter('SingleCSAT') Call EigenSystem(ux, ux, vx, vx, Rinv, Rdir, lamda) if ( cleslog_active(CLESLOG_CBC, CLESLOG_CBC_EIGENSYSTEM) & .eq. CLES_TRUE ) then call cles_test_pack_dir(ipos, jpos, kpos) call cles_test_eigen(ux, Rinv, Rdir, lamda) call cleslog_log_flush() endif c = (lamda(nvars)-lamda(1))*0.5d0 ! for example DeltaG(:) = 0.0 do k=1,3 do i=1,nvars vec(i) = SUM(Rinv(i,:)*flux(:,k)) enddo flux(:,k) = vec(:) enddo if ( useSource .eq. CLES_TRUE ) then ! do i=1, nvars ! dU(i) = SUM(Rinv(i,:)*source(:)) ! enddo ! source = dU endif ! save characteristic incomming strength fx = flux(1,1) fy = flux(1,2) fz = flux(1,3) flux(:,:) = 0.0d0 ! Point #1 ! The only values that should be non-zero are those values that ! correspond to characteristic waves trying to get inside the domain ! for which we must know their value. if ( lamda(nvars) .lt. 0.0 ) then DeltaG(:) = lamda(:) ! if ( useSource .eq. CLES_TRUE ) source(:) = 0.0d0 else if ( lamda(2) .lt. 0.0 ) then DeltaG(1:nvars-1) = lamda(1:nvars-1) ! if ( useSource .eq. CLES_TRUE ) source(1:nvars-1) = 0.0d0 else if ( lamda(1) .lt. 0.0 ) then DeltaG(1) = lamda(1) ! if ( useSource .eq. CLES_TRUE ) source(1) = 0.0d0 else ! Nothing to do endif a = tcd_bndry(1,1)/dl ! Add penalty terms for those directions that are incomming dU(1:nvars) = ux(1:nvars)-uxref(1:nvars) do i=1, nvars vec(i) = DeltaG(i)*a*tau(i)*SUM(Rinv(i,:)*dU(:)) enddo if ( cbc_type .eq. CLES_CBC_OUTFLOW .and. lamda(1) .lt. 0.0d0 & .and. cbc_Mach .lt. 1.0d0 ) then if ( cleslog_active(CLESLOG_CBC, CLESLOG_CBC_REFPRESSURE) & .eq. CLES_TRUE ) then call cles_test_pack_dir(ipos, jpos, kpos) write(slog, '(A10,E12.4)') 'Pressure =', p call cleslog_log(slog//char(0)) call cleslog_log_flush() endif vec(1) = 0.0d0 ! if ( useSource .eq. CLES_TRUE ) source(1) = 0.0d0 ! cancel outgoing differential part flux(1,1) = cbc_factor*(vx(5)-p)/((xrg-xlg)*c*Rdir(1,1))-fx if ( cbc_mode .eq. CLES_CBC_MODE2D ) then flux(1,2) = - fy flux(1,3) = - fz endif endif if ( useSource .eq. CLES_TRUE ) then ! do i=1, nvars ! dU(i) = SUM(Rdir(i,:)*source(:)) ! enddo ! source = dU endif ! Multiply by R (make conservative) do i=1, nvars sat(i) = SUM(Rdir(i,:)*vec(:)) enddo ! rotate back flux part do k=1,3 do i=1,nvars vec(i) = SUM(Rdir(i,:)*flux(:,k)) enddo flux(:,k) = vec(:) enddo call cleslog_log_exit('SingleCSAT') END SUBROUTINE SingleCSAT subroutine DirectRotateXYZ(ux, vx, uxref, src, flux, P, bc) use method_parms use Generic_Source implicit none include 'cles.i' double precision, intent(inout) :: ux(ncomps), vx(nvars) double precision, intent(inout) :: uxref(ncomps), src(nvars) double precision, intent(inout) :: flux(nvars,3) double precision, intent(in) :: P(3,3) integer, intent(in) :: bc integer i, k double precision vec(3) ! vector rotation do i=1,3 vec(i) = SUM(ux(2:4)*P(:,i)) enddo ux(2:4) = vec(:) do i=1,3 vec(i) = SUM(vx(2:4)*P(:,i)) enddo vx(2:4) = vec(:) do i=1,3 vec(i) = SUM(uxref(2:4)*P(:,i)) enddo uxref(2:4) = vec(:) if ( useSource .eq. CLES_TRUE ) then do i=1,3 vec(i) = SUM(src(2:4)*P(:,i)) enddo src(2:4) = vec(:) endif ! change of independent variables do i=1,nvars do k=1,3 vec(k) = SUM(P(:,k)*flux(i,:)) enddo flux(i,:) = vec(:) enddo ! rotation do k=1, 3 do i=1,3 vec(i) = SUM(flux(2:4,k)*P(:,i)) enddo flux(2:4,k) = vec(:) enddo end subroutine DirectRotateXYZ subroutine InverseRotateXYZ(src, sat, flux, P, bc) use method_parms use Generic_Source implicit none include 'cles.i' double precision, intent(inout) :: src(nvars), sat(nvars) double precision, intent(inout) :: flux(nvars,3) double precision, intent(in) :: P(3,3) integer :: bc double precision vec(3) integer i, k do i=1,3 vec(i) = SUM(P(i,:)*sat(2:4)) enddo sat(2:4) = vec(:) if ( useSource .eq. CLES_TRUE ) then do i=1,3 vec(i) = SUM(P(i,:)*src(2:4)) enddo src(2:4) = vec(:) endif ! rotation first do k=1, 3 do i=1,3 vec(i) = SUM(P(i,:)*flux(2:4,k)) enddo flux(2:4,k) = vec(:) enddo ! change of independent variables do i=1,nvars do k=1,3 vec(k) = SUM(flux(i,:)*P(k,:)) enddo flux(i,:) = vec(:) enddo end subroutine InverseRotateXYZ subroutine setproj(theta, phi, proj) implicit none double precision theta, phi, proj(3,3) double precision ctheta, stheta, cphi, sphi ! projection matrix ctheta = cos(theta) stheta = sin(theta) cphi = cos(phi) sphi = sin(phi) proj(1,1) = ctheta*cphi proj(1,2) = -stheta proj(1,3) = -sphi*ctheta proj(2,1) = cphi*stheta proj(2,2) = ctheta proj(2,3) = -sphi*stheta proj(3,1) = sphi proj(3,2) = 0.0d0 proj(3,3) = cphi return end subroutine setproj SUBROUTINE OneDAcousticBC(fxi, ux, vx, source, dcflag) ! ---- Shared Variables USE mesh USE array_bounds USE method_parms USE Generic_Source IMPLICIT NONE include 'cles.i' DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi) DOUBLE PRECISION,INTENT(INOUT) :: fxi(ncomps,ixlo:ixhi,1) ! the cell centered flux DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi) DOUBLE PRECISION, INTENT(INOUT) :: source(nvars,1:nx) INTEGER, INTENT(IN) :: dcflag(1:nx+1,1) DOUBLE PRECISION :: uxloc(ncomps), vxloc(nvars), uxref(ncomps) DOUBLE PRECISION :: sat(nvars), flux(nvars,3), Proj(3,3) if ( cbc_active .ne. CLES_TRUE ) return !----------------------------------------------------------- ! ! LEFT BOUNDARY ! !----------------------------------------------------------- if ( cbc_ixlow .ne. CLES_CBC_NONE ) then Proj = 0.0d0 Proj(1,1) = -1.0d0 Proj(2,2) = -1.0d0 Proj(3,3) = 1.0d0 flux(:,1) = (fxi(1:nvars,2,1)-fxi(1:nvars,1,1))/dx flux(:,2:3) = 0.0d0 uxloc = ux(:,1) vxloc = vx(:,1) uxref = ux(:,0) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(1,1), flux, Proj, & cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(1,1), sat, flux, vx(5,0), & dx, cbc_force_matrix(:,1), cbc_ixlow, 1, 0, 0) ! rotation back call InverseRotateXYZ(source(1,1), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1) = source(:,1) - sat(:) - flux(:,1) endif !----------------------------------------------------------- ! ! RIGHT BOUNDARY ! !----------------------------------------------------------- if ( cbc_ixup .ne. CLES_CBC_NONE ) then flux(:,1) = (fxi(1:nvars,nx+1,1)-fxi(1:nvars,nx,1))/dx flux(:,2:3) = 0.0d0 call SingleCSAT(ux(1,nx), vx(1,nx), ux(1,nx+1), source(1,nx), sat, & flux, vx(5,nx+1), dx, cbc_force_matrix(:,2), cbc_ixup, nx, 0, 0) ! Correct flux to the last point source(:,nx) = source(:,nx) - sat(:) - flux(:,1) endif RETURN END SUBROUTINE OneDAcousticBC SUBROUTINE TwoDAcousticBC(fxi, ux, vx, source, dcflag) ! ---- Shared Variables USE mesh USE array_bounds USE method_parms USE Generic_Source IMPLICIT NONE include 'cles.i' DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi) DOUBLE PRECISION, INTENT(INOUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,2) ! the cell centered flux DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi) DOUBLE PRECISION, INTENT(INOUT) :: source(nvars,1:nx,1:ny) INTEGER, INTENT(IN) :: dcflag(1:nx+1,1:ny+1,2) INTEGER :: i, j DOUBLE PRECISION :: uxloc(ncomps), vxloc(nvars), uxref(ncomps) DOUBLE PRECISION :: Proj(3,3), flux(nvars,3), sat(nvars) INTEGER :: ilow, iup, jlow, jup DOUBLE PRECISION :: p, dl, isqr2 if ( cbc_active .ne. CLES_TRUE ) return isqr2 = 1.0d0/sqrt(2.0d0) !----------------------------------------------------------- ! ! LEFT BOUNDARY ! !----------------------------------------------------------- if ( cbc_ixlow .ne. CLES_CBC_NONE ) then jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_ixlow ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_ixlow ) jup = 1 ! low priority corner (just avoid) if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 endif Proj = 0.0d0 Proj(1,1) = -1.0d0 Proj(2,2) = -1.0d0 Proj(3,3) = 1.0d0 do j=1+jlow, ny-jup flux(:,1) = (fxi(1:nvars,2,j,1) -fxi(1:nvars,1,j,1))/dx flux(:,2) = (fxi(1:nvars,1,j+1,2)-fxi(1:nvars,1,j,2))/dy flux(:,3) = 0.0d0 uxloc = ux(:,1,j) vxloc = vx(:,1,j) uxref = ux(:,0,j) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,j), & flux, Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,j), sat, & flux, vx(5,0,j), dx, cbc_force_matrix(:,1), cbc_ixlow, & 1, j, 0) ! rotation back call InverseRotateXYZ(source(:,1,j), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,j) = source(:,1,j) - sat(:) - flux(:,1) - flux(:,2) enddo endif !----------------------------------------------------------- ! ! RIGHT BOUNDARY ! !----------------------------------------------------------- if ( cbc_ixup .ne. CLES_CBC_NONE ) then jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_ixup ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_ixup ) jup = 1 ! low priority corner (just avoid) if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 endif do j=1+jlow, ny-jup flux(:,1) = (fxi(1:nvars,nx+1,j,1)-fxi(1:nvars,nx,j,1))/dx flux(:,2) = (fxi(1:nvars,nx,j+1,2)-fxi(1:nvars,nx,j,2))/dy flux(:,3) = 0.0d0 call SingleCSAT(ux(:,nx,j), vx(:,nx,j), ux(:,nx+1,j), & source(:,nx,j), sat, flux, vx(5,nx+1,j), dx, & cbc_force_matrix(:,2), cbc_ixup, nx, j, 0) ! Correct flux to the last point source(:,nx,j) = source(:,nx,j) - sat(:) -flux(:,1) -flux(:,2) enddo endif !----------------------------------------------------------- ! ! BOTTOM BOUNDARY ! !----------------------------------------------------------- if ( cbc_iylow .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_iylow ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_iylow ) iup = 1 ! low priority corner (just avoid) if ( cbc_iylow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 endif Proj = 0.0d0 Proj(1,2) = 1.0d0 Proj(2,1) = -1.0d0 Proj(3,3) = 1.0d0 do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,1,1)-fxi(1:nvars,i,1,1))/dx flux(:,2) = (fxi(1:nvars,i,2,2) -fxi(1:nvars,i,1,2))/dy flux(:,3) = 0.0d0 uxloc = ux(:,i,1) vxloc = vx(:,i,1) uxref = ux(:,i,0) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,1), flux, & Proj, cbc_iylow) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,1), sat, & flux, vx(5,i,0), dy, cbc_force_matrix(:,3), cbc_iylow, & i, 1, 0) ! rotation back call InverseRotateXYZ(source(:,i,1), sat, flux, Proj, cbc_iylow) ! Correct flux to the last point source(:,i,1) = source(:,i,1) - sat(:)- flux(:,1)- flux(:,2) enddo endif !----------------------------------------------------------- ! ! TOP BOUNDARY ! !----------------------------------------------------------- if ( cbc_iyup .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_iyup ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_iyup ) iup = 1 ! low priority corner (just avoid) if ( cbc_iyup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 endif Proj = 0.0d0 Proj(1,2) = -1.0d0 Proj(2,1) = 1.0d0 Proj(3,3) = 1.0d0 do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,ny,1)-fxi(1:nvars,i,ny,1))/dx flux(:,2) = (fxi(1:nvars,i,ny+1,2)-fxi(1:nvars,i,ny,2))/dy flux(:,3) = 0.0d0 uxloc = ux(:,i,ny) vxloc = vx(:,i,ny) uxref = ux(:,i,ny+1) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,ny), & flux, Proj, cbc_iyup) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,ny), sat, & flux, vx(5,i,ny+1), dy, cbc_force_matrix(:,4), & cbc_iyup, i, ny, 0) ! rotation back call InverseRotateXYZ(source(:,i,ny), sat, flux, Proj, cbc_iyup) ! Correct flux to the last point source(:,i,ny) = source(:,i,ny) - sat(:)- flux(:,1)- flux(:,2) enddo endif !----------------------------------------------------------- ! ! LOWER LEFT CORNER ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_iylow .and. cbc_ixlow .ne. CLES_CBC_NONE ) then Proj = 0.0d0 Proj(1,1) = -isqr2 Proj(1,2) = isqr2 Proj(2,1) = -isqr2 Proj(2,2) = -isqr2 Proj(3,3) = 1.0d0 flux(:,1) = (fxi(1:nvars,2,1,1)-fxi(1:nvars,1,1,1))/dx flux(:,2) = (fxi(1:nvars,1,2,2)-fxi(1:nvars,1,1,2))/dy flux(:,3) = 0.0d0 uxloc = ux(:,1,1) vxloc = vx(:,1,1) uxref = (ux(:,1,0)+ux(:,0,1))*0.5d0 dl = sqrt(dx*dy) p = (vx(5,1,0)+vx(5,0,1))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,1), flux, & Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,1), sat, flux, & p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, 1, 0) ! rotation back call InverseRotateXYZ(source(:,1,1), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,1) = source(:,1,1) - sat(:)- flux(:,1) -flux(:,2) endif !----------------------------------------------------------- ! ! UPPER LEFT CORNER ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_iyup .and. cbc_ixlow .ne. CLES_CBC_NONE ) then Proj = 0.0d0 Proj(1,1) = -isqr2 Proj(1,2) = -isqr2 Proj(2,1) = isqr2 Proj(2,2) = -isqr2 Proj(3,3) = 1.0d0 flux(:,1) = (fxi(1:nvars,2,ny,1) -fxi(1:nvars,1,ny,1))/dx flux(:,2) = (fxi(1:nvars,1,ny+1,2)-fxi(1:nvars,1,ny,2))/dy flux(:,3) = 0.0d0 uxloc = ux(:,1,ny) vxloc = vx(:,1,ny) uxref = (ux(:,1,ny+1)+ux(:,0,ny))*0.5d0 dl = sqrt(dx*dy) p = (vx(5,1,ny+1)+vx(5,0,ny))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,ny), flux, & Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,ny), sat, flux, & p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, ny, 0) ! rotation back call InverseRotateXYZ(source(:,1,ny), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,ny) = source(:,1,ny) - sat(:)- flux(:,1) -flux(:,2) endif !----------------------------------------------------------- ! ! LOWER RIGHT CORNER ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_iylow .and. cbc_ixup .ne. CLES_CBC_NONE ) then Proj = 0.0d0 Proj(1,1) = isqr2 Proj(1,2) = isqr2 Proj(2,1) = -isqr2 Proj(2,2) = isqr2 Proj(3,3) = 1.0d0 flux(:,1) = (fxi(1:nvars,nx+1,1,1)-fxi(1:nvars,nx,1,1))/dx flux(:,2) = (fxi(1:nvars,nx,2,2) -fxi(1:nvars,nx,1,2))/dy flux(:,3) = 0.0d0 uxloc = ux(:,nx,1) vxloc = vx(:,nx,1) uxref = (ux(:,nx,0)+ux(:,nx+1,1))*0.5d0 dl = sqrt(dx*dy) p = (vx(5,nx,0)+vx(5,nx+1,1))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,1), flux, & Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,1), sat, flux, & p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, 1, 0) ! rotation back call InverseRotateXYZ(source(:,nx,1), sat, flux, Proj, cbc_ixup) ! Correct flux to the last point source(:,nx,1) = source(:,nx,1) - sat(:)- flux(:,1) -flux(:,2) endif !----------------------------------------------------------- ! ! UPPER RIGHT CORNER ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_iyup .and. cbc_ixup .ne. CLES_CBC_NONE ) then Proj = 0.0d0 Proj(1,1) = isqr2 Proj(1,2) = -isqr2 Proj(2,1) = isqr2 Proj(2,2) = isqr2 Proj(3,3) = 1.0d0 flux(:,1) = (fxi(1:nvars,nx+1,ny,1)-fxi(1:nvars,nx,ny,1))/dx flux(:,2) = (fxi(1:nvars,nx,ny+1,2)-fxi(1:nvars,nx,ny,2))/dy flux(:,3) = 0.0d0 uxloc = ux(:,nx,ny) vxloc = vx(:,nx,ny) uxref = (ux(:,nx,ny+1)+ux(:,nx+1,ny))*0.5d0 dl = sqrt(dx*dy) p = (vx(5,nx,ny+1)+vx(5,nx+1,ny))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,ny), flux, & Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,ny), sat, flux, & p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, ny, 0) ! rotation back call InverseRotateXYZ(source(:,nx,ny), sat, flux, Proj, cbc_ixup) ! Correct flux to the last point source(:,nx,ny) = source(:,nx,ny) - sat(:)- flux(:,1) - flux(:,2) endif RETURN END SUBROUTINE TwoDAcousticBC SUBROUTINE ThreeDAcousticBC(fxi, ux, vx, source, dcflag) ! ---- Shared Variables USE mesh USE array_bounds USE method_parms USE Generic_Source IMPLICIT NONE include 'cles.i' DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi) DOUBLE PRECISION, INTENT(INOUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi,3) ! the cell centered flux DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi) DOUBLE PRECISION, INTENT(INOUT) :: source(nvars,1:nx,1:ny,1:nz) INTEGER, INTENT(IN) :: dcflag(1:nx+1,1:ny+1,1:nz+1,3) INTEGER :: i, j, k INTEGER :: ilow, iup, jlow, jup, klow, kup DOUBLE PRECISION :: uxloc(ncomps), vxloc(nvars), uxref(ncomps) DOUBLE PRECISION :: sat(nvars), flux(nvars,3), Proj(3,3), dl, p DOUBLE PRECISION :: isqr2, pi, phi if ( cbc_active .ne. CLES_TRUE ) return isqr2 = 1.0d0/sqrt(2.0d0) pi = 4.0d0*atan(1.0d0) phi = asin(1.0d0/sqrt(3.0d0)) !----------------------------------------------------------- ! ! LEFT BOUNDARY ! !----------------------------------------------------------- if ( cbc_ixlow .ne. CLES_CBC_NONE ) then jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_ixlow ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_ixlow ) jup = 1 klow = 0; kup = 0 ! lower corner if ( cbc_izlow .eq. cbc_ixlow ) klow = 1 ! upper corner if ( cbc_izup .eq. cbc_ixlow ) kup = 1 ! low priority corner (just avoid) if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1 if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1 endif Proj = 0.0d0 Proj(1,1) = -1.0d0 Proj(2,2) = -1.0d0 Proj(3,3) = 1.0d0 do k=1+klow, nz-kup do j=1+jlow, ny-jup flux(:,1) = (fxi(1:nvars,2,j,k,1) -fxi(1:nvars,1,j,k,1))/dx flux(:,2) = (fxi(1:nvars,1,j+1,k,2)-fxi(1:nvars,1,j,k,2))/dy flux(:,3) = (fxi(1:nvars,1,j,k+1,3)-fxi(1:nvars,1,j,k,3))/dz uxloc = ux(:,1,j,k) vxloc = vx(:,1,j,k) uxref = ux(:,0,j,k) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,j,k), & flux, Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,j,k), sat, & flux, vx(5,0,j,k), dx, cbc_force_matrix(:,1), cbc_ixlow, & 1, j, k) ! rotation back call InverseRotateXYZ(source(:,1,j,k), sat, flux, Proj, & cbc_ixlow) ! Correct flux to the last point source(:,1,j,k) = source(:,1,j,k) - sat(:) & -flux(:,1) -flux(:,2) -flux(:,3) enddo enddo endif !----------------------------------------------------------- ! ! RIGHT BOUNDARY ! !----------------------------------------------------------- if ( cbc_ixup .ne. CLES_CBC_NONE ) then jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_ixup ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_ixup ) jup = 1 klow = 0; kup = 0 ! lower corner if ( cbc_izlow .eq. cbc_ixup ) klow = 1 ! upper corner if ( cbc_izup .eq. cbc_ixup ) kup = 1 ! low priority corner (just avoid) if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1 if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1 endif do k=1+klow, nz-kup do j=1+jlow, ny-jup flux(:,1) = (fxi(1:nvars,nx+1,j,k,1)-fxi(1:nvars,nx,j,k,1))/dx flux(:,2) = (fxi(1:nvars,nx,j+1,k,2)-fxi(1:nvars,nx,j,k,2))/dy flux(:,3) = (fxi(1:nvars,nx,j,k+1,3)-fxi(1:nvars,nx,j,k,3))/dz call SingleCSAT(ux(:,nx,j,k), vx(:,nx,j,k), ux(:,nx+1,j,k), & source(:,nx,j,k), sat, flux, vx(5,nx+1,j,k), dx, & cbc_force_matrix(:,2), cbc_ixup, nx, j, k) ! Correct flux to the last point source(:,nx,j,k) = source(:,nx,j,k) - sat(:) & - flux(:,1)- flux(:,2)- flux(:,3) enddo enddo endif !----------------------------------------------------------- ! ! BOTTOM BOUNDARY ! !----------------------------------------------------------- if ( cbc_iylow .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_iylow ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_iylow ) iup = 1 klow = 0; kup = 0 ! lower corner if ( cbc_izlow .eq. cbc_iylow ) klow = 1 ! upper corner if ( cbc_izup .eq. cbc_iylow ) kup = 1 ! low priority corner (just avoid) if ( cbc_iylow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1 if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1 endif Proj = 0.0d0 Proj(1,2) = 1.0d0 Proj(2,1) = -1.0d0 Proj(3,3) = 1.0d0 do k=1+klow, nz-kup do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,1,k,1)-fxi(1:nvars,i,1,k,1))/dx flux(:,2) = (fxi(1:nvars,i,2,k,2) -fxi(1:nvars,i,1,k,2))/dy flux(:,3) = (fxi(1:nvars,i,1,k+1,3)-fxi(1:nvars,i,1,k,3))/dz uxloc = ux(:,i,1,k) vxloc = vx(:,i,1,k) uxref = ux(:,i,0,k) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,1,k), & flux, Proj, cbc_iylow) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,1,k), sat, & flux, vx(5,i,0,k), dy, cbc_force_matrix(:,3), cbc_iylow, & i, 1, k) ! rotation back call InverseRotateXYZ(source(:,i,1,k), sat, flux, Proj, & cbc_iylow) ! Correct flux to the last point source(:,i,1,k) = source(:,i,1,k) - sat(:)- flux(:,1)& - flux(:,2)- flux(:,3) enddo enddo endif !----------------------------------------------------------- ! ! TOP BOUNDARY ! !----------------------------------------------------------- if ( cbc_iyup .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_iyup ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_iyup ) iup = 1 klow = 0; kup = 0 ! lower corner if ( cbc_izlow .eq. cbc_iyup ) klow = 1 ! upper corner if ( cbc_izup .eq. cbc_iyup ) kup = 1 ! low priority corner (just avoid) if ( cbc_iyup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1 if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1 endif Proj = 0.0d0 Proj(1,2) = -1.0d0 Proj(2,1) = 1.0d0 Proj(3,3) = 1.0d0 do k=1+klow, nz-kup do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,ny,k,1)-fxi(1:nvars,i,ny,k,1))/dx flux(:,2) = (fxi(1:nvars,i,ny+1,k,2)-fxi(1:nvars,i,ny,k,2))/dy flux(:,3) = (fxi(1:nvars,i,ny,k+1,3)-fxi(1:nvars,i,ny,k,3))/dz uxloc = ux(:,i,ny,k) vxloc = vx(:,i,ny,k) uxref = ux(:,i,ny+1,k) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,ny,k), & flux, Proj, cbc_iyup) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,ny,k), sat, & flux, vx(5,i,ny+1,k), dy, cbc_force_matrix(:,4), & cbc_iyup, i, ny, k) ! rotation back call InverseRotateXYZ(source(:,i,ny,k), sat, flux, Proj, & cbc_iyup) ! Correct flux to the last point source(:,i,ny,k) = source(:,i,ny,k) - sat(:)- flux(:,1)& - flux(:,2)- flux(:,3) enddo enddo endif !----------------------------------------------------------- ! ! BACK BOUNDARY ! !----------------------------------------------------------- if ( cbc_izlow .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_izlow ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_izlow ) iup = 1 jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_izlow ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_izlow ) jup = 1 ! low priority corner (just avoid) if ( cbc_izlow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 endif Proj = 0.0d0 Proj(1,3) = 1.0d0 Proj(2,2) = 1.0d0 Proj(3,1) = -1.0d0 do j=1+jlow, ny-jup do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,j,1,1)-fxi(1:nvars,i,j,1,1))/dx flux(:,2) = (fxi(1:nvars,i,j+1,1,2)-fxi(1:nvars,i,j,1,2))/dy flux(:,3) = (fxi(1:nvars,i,j,2,3) -fxi(1:nvars,i,j,1,3))/dz uxloc = ux(:,i,j,1) vxloc = vx(:,i,j,1) uxref = ux(:,i,j,0) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,j,1), & flux, Proj, cbc_izlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,j,1), sat, & flux, vx(5,i,j,0), dz, cbc_force_matrix(:,5), cbc_izlow, & i, j, 1) ! rotation back call InverseRotateXYZ(source(:,i,j,1), sat, flux, Proj, & cbc_izlow) ! Correct flux to the last point source(:,i,j,1) = source(:,i,j,1) - sat(:)- flux(:,1)& - flux(:,2)- flux(:,3) enddo enddo endif !----------------------------------------------------------- ! ! FRONT BOUNDARY ! !----------------------------------------------------------- if ( cbc_izup .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_izup ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_izup ) iup = 1 jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_izup ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_izup ) jup = 1 ! low priority corner (just avoid) if ( cbc_izup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 endif Proj = 0.0d0 Proj(1,3) = -1.0d0 Proj(2,2) = 1.0d0 Proj(3,1) = 1.0d0 do j=1+jlow, ny-jup do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,j,nz,1)-fxi(1:nvars,i,j,nz,1))/dx flux(:,2) = (fxi(1:nvars,i,j+1,nz,2)-fxi(1:nvars,i,j,nz,2))/dy flux(:,3) = (fxi(1:nvars,i,j,nz+1,3)-fxi(1:nvars,i,j,nz,3))/dz uxloc = ux(:,i,j,nz) vxloc = vx(:,i,j,nz) uxref = ux(:,i,j,nz+1) ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,j,nz), & flux, Proj, cbc_izup) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,j,nz), sat, & flux, vx(5,i,j,nz+1), dz, cbc_force_matrix(:,6), & cbc_izup, i, j, nz) ! rotation back call InverseRotateXYZ(source(:,i,j,nz), sat, flux, Proj, & cbc_izup) ! Correct flux to the last point source(:,i,j,nz) = source(:,i,j,nz) - sat(:)- flux(:,1)& - flux(:,2)- flux(:,3) enddo enddo endif !----------------------------------------------------------- ! ! LOWER LEFT EDGE ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_iylow .and. cbc_ixlow .ne. CLES_CBC_NONE ) then klow = 0; kup = 0 ! lower corner if ( cbc_izlow .eq. cbc_ixlow ) klow = 1 ! upper corner if ( cbc_izup .eq. cbc_ixlow ) kup = 1 ! low priority corner (just avoid) if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1 if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1 endif Proj = 0.0d0 Proj(1,1) = -isqr2 Proj(1,2) = isqr2 Proj(2,1) = -isqr2 Proj(2,2) = -isqr2 Proj(3,3) = 1.0d0 dl = sqrt(dx*dy) do k=1+klow, nz-kup flux(:,1) = (fxi(1:nvars,2,1,k,1)-fxi(1:nvars,1,1,k,1))/dx flux(:,2) = (fxi(1:nvars,1,2,k,2)-fxi(1:nvars,1,1,k,2))/dy flux(:,3) = (fxi(1:nvars,1,1,k+1,3)-fxi(1:nvars,1,1,k,3))/dz uxloc = ux(:,1,1,k) vxloc = vx(:,1,1,k) uxref = (ux(:,1,0,k)+ux(:,0,1,k))*0.5d0 p = (vx(5,1,0,k)+vx(5,0,1,k))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,1,k), flux, & Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,1,k), sat, & flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, 1, k) ! rotation back call InverseRotateXYZ(source(:,1,1,k), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,1,k) = source(:,1,1,k) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) enddo endif !----------------------------------------------------------- ! ! UPPER LEFT EDGE ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_iyup .and. cbc_ixlow .ne. CLES_CBC_NONE ) then klow = 0; kup = 0 ! lower corner if ( cbc_izlow .eq. cbc_ixlow ) klow = 1 ! upper corner if ( cbc_izup .eq. cbc_ixlow ) kup = 1 ! low priority corner (just avoid) if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1 if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1 endif Proj = 0.0d0 Proj(1,1) = -isqr2 Proj(1,2) = -isqr2 Proj(2,1) = isqr2 Proj(2,2) = -isqr2 Proj(3,3) = 1.0d0 dl = sqrt(dx*dy) do k=1+klow, nz-kup flux(:,1) = (fxi(1:nvars,2,ny,k,1) -fxi(1:nvars,1,ny,k,1))/dx flux(:,2) = (fxi(1:nvars,1,ny+1,k,2)-fxi(1:nvars,1,ny,k,2))/dy flux(:,3) = (fxi(1:nvars,1,ny,k+1,3)-fxi(1:nvars,1,ny,k,3))/dz uxloc = ux(:,1,ny,k) vxloc = vx(:,1,ny,k) uxref = (ux(:,1,ny+1,k)+ux(:,0,ny,k))*0.5d0 p = (vx(5,1,ny+1,k)+vx(5,0,ny,k))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,ny,k), flux,& Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,ny,k), sat, & flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, ny, k) ! rotation back call InverseRotateXYZ(source(:,1,ny,k), sat, flux, Proj, & cbc_ixlow) ! Correct flux to the last point source(:,1,ny,k) = source(:,1,ny,k) - sat(:) & - flux(:,1) -flux(:,2) -flux(:,3) enddo endif !----------------------------------------------------------- ! ! LOWER RIGHT EDGE ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_iylow .and. cbc_ixup .ne. CLES_CBC_NONE ) then klow = 0; kup = 0 ! lower corner if ( cbc_izlow .eq. cbc_ixup ) klow = 1 ! upper corner if ( cbc_izup .eq. cbc_ixup ) kup = 1 ! low priority corner (just avoid) if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1 if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1 endif Proj = 0.0d0 Proj(1,1) = isqr2 Proj(1,2) = isqr2 Proj(2,1) = -isqr2 Proj(2,2) = isqr2 Proj(3,3) = 1.0d0 dl = sqrt(dx*dy) do k=1+klow, nz-kup flux(:,1) = (fxi(1:nvars,nx+1,1,k,1)-fxi(1:nvars,nx,1,k,1))/dx flux(:,2) = (fxi(1:nvars,nx,2,k,2) -fxi(1:nvars,nx,1,k,2))/dy flux(:,3) = (fxi(1:nvars,nx,1,k+1,3)-fxi(1:nvars,nx,1,k,3))/dz uxloc = ux(:,nx,1,k) vxloc = vx(:,nx,1,k) uxref = (ux(:,nx,0,k)+ux(:,nx+1,1,k))*0.5d0 p = (vx(5,nx,0,k)+vx(5,nx+1,1,k))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,1,k), flux,& Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,1,k), sat, & flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, 1, k) ! rotation back call InverseRotateXYZ(source(:,nx,1,k), sat, flux, Proj, cbc_ixup) ! Correct flux to the last point source(:,nx,1,k) = source(:,nx,1,k) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) enddo endif !----------------------------------------------------------- ! ! UPPER RIGHT EDGE ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_iyup .and. cbc_ixup .ne. CLES_CBC_NONE ) then klow = 0; kup = 0 ! lower corner if ( cbc_izlow .eq. cbc_ixup ) klow = 1 ! upper corner if ( cbc_izup .eq. cbc_ixup ) kup = 1 ! low priority corner (just avoid) if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1 if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1 endif Proj = 0.0d0 Proj(1,1) = isqr2 Proj(1,2) = -isqr2 Proj(2,1) = isqr2 Proj(2,2) = isqr2 Proj(3,3) = 1.0d0 dl = sqrt(dx*dy) do k=1+klow, nz-kup flux(:,1) = (fxi(1:nvars,nx+1,ny,k,1)-fxi(1:nvars,nx,ny,k,1))/dx flux(:,2) = (fxi(1:nvars,nx,ny+1,k,2)-fxi(1:nvars,nx,ny,k,2))/dy flux(:,3) = (fxi(1:nvars,nx,ny,k+1,3)-fxi(1:nvars,nx,ny,k,3))/dz uxloc = ux(:,nx,ny,k) vxloc = vx(:,nx,ny,k) uxref = (ux(:,nx,ny+1,k)+ux(:,nx+1,ny,k))*0.5d0 p = (vx(5,nx,ny+1,k)+vx(5,nx+1,ny,k))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,ny,k), & flux, Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,ny,k), sat, & flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, ny, k) ! rotation back call InverseRotateXYZ(source(:,nx,ny,k), sat, flux, Proj, & cbc_ixup) ! Correct flux to the last point source(:,nx,ny,k) = source(:,nx,ny,k) - sat(:) & - flux(:,1) - flux(:,2) -flux(:,3) enddo endif !----------------------------------------------------------- ! ! BACK LEFT EDGE ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_izlow .and. cbc_ixlow .ne. CLES_CBC_NONE ) then jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_ixlow ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_ixlow ) jup = 1 ! low priority corner (just avoid) if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 endif Proj = 0.0d0 Proj(1,1) = -isqr2 Proj(1,3) = isqr2 Proj(3,1) = -isqr2 Proj(3,3) = -isqr2 Proj(2,2) = 1.0d0 dl = sqrt(dx*dz) do j=1+jlow, ny-jup flux(:,1) = (fxi(1:nvars,2,j,1,1) -fxi(1:nvars,1,j,1,1))/dx flux(:,2) = (fxi(1:nvars,1,j+1,1,2)-fxi(1:nvars,1,j,1,2))/dy flux(:,3) = (fxi(1:nvars,1,j,2,3) -fxi(1:nvars,1,j,1,3))/dz uxloc = ux(:,1,j,1) vxloc = vx(:,1,j,1) uxref = (ux(:,1,j,0)+ux(:,0,j,1))*0.5d0 p = (vx(5,1,j,0)+vx(5,0,j,1))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,j,1), flux, & Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,j,1), sat, & flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, j, 1) ! rotation back call InverseRotateXYZ(source(:,1,j,1), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,j,1) = source(:,1,j,1) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) enddo endif !----------------------------------------------------------- ! ! FRONT LEFT EDGE ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_izup .and. cbc_ixlow .ne. CLES_CBC_NONE ) then jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_ixlow ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_ixlow ) jup = 1 ! low priority corner (just avoid) if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 endif Proj = 0.0d0 Proj(1,1) = -isqr2 Proj(1,3) = -isqr2 Proj(3,1) = isqr2 Proj(3,3) = -isqr2 Proj(2,2) = 1.0d0 dl = sqrt(dx*dz) do j=1+jlow, ny-jup flux(:,1) = (fxi(1:nvars,2,j,nz,1) -fxi(1:nvars,1,j,nz,1))/dx flux(:,2) = (fxi(1:nvars,1,j+1,nz,2)-fxi(1:nvars,1,j,nz,2))/dy flux(:,3) = (fxi(1:nvars,1,j,nz+1,3)-fxi(1:nvars,1,j,nz,3))/dz uxloc = ux(:,1,j,nz) vxloc = vx(:,1,j,nz) uxref = (ux(:,1,j,nz+1)+ux(:,0,j,nz))*0.5d0 p = (vx(5,1,j,nz+1)+vx(5,0,j,nz))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,j,nz), flux,& Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,j,nz), sat, & flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, j, nz) ! rotation back call InverseRotateXYZ(source(:,1,j,nz), sat, flux, Proj, & cbc_ixlow) ! Correct flux to the last point source(:,1,j,nz) = source(:,1,j,nz) - sat(:) & - flux(:,1) -flux(:,2) -flux(:,3) enddo endif !----------------------------------------------------------- ! ! BACK RIGHT EDGE ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_izlow .and. cbc_ixup .ne. CLES_CBC_NONE ) then jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_ixup ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_ixup ) jup = 1 ! low priority corner (just avoid) if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 endif Proj = 0.0d0 Proj(1,1) = isqr2 Proj(1,3) = isqr2 Proj(3,1) = -isqr2 Proj(3,3) = isqr2 Proj(2,2) = 1.0d0 dl = sqrt(dx*dz) do j=1+jlow, ny-jup flux(:,1) = (fxi(1:nvars,nx+1,j,1,1)-fxi(1:nvars,nx,j,1,1))/dx flux(:,2) = (fxi(1:nvars,nx,j+1,1,2)-fxi(1:nvars,nx,j,1,2))/dy flux(:,3) = (fxi(1:nvars,nx,j,2,3) -fxi(1:nvars,nx,j,1,3))/dz uxloc = ux(:,nx,j,1) vxloc = vx(:,nx,j,1) uxref = (ux(:,nx,j,0)+ux(:,nx+1,j,1))*0.5d0 p = (vx(5,nx,j,0)+vx(5,nx+1,j,1))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,j,1), flux,& Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,j,1), sat, & flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, j, 1) ! rotation back call InverseRotateXYZ(source(:,nx,j,1), sat, flux, Proj, cbc_ixup) ! Correct flux to the last point source(:,nx,j,1) = source(:,nx,j,1) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) enddo endif !----------------------------------------------------------- ! ! FRONT RIGHT EDGE ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_izup .and. cbc_ixup .ne. CLES_CBC_NONE ) then jlow = 0; jup = 0 ! lower corner if ( cbc_iylow .eq. cbc_ixup ) jlow = 1 ! upper corner if ( cbc_iyup .eq. cbc_ixup ) jup = 1 ! low priority corner (just avoid) if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1 if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1 endif Proj = 0.0d0 Proj(1,1) = isqr2 Proj(1,3) = -isqr2 Proj(3,1) = isqr2 Proj(3,3) = isqr2 Proj(2,2) = 1.0d0 dl = sqrt(dx*dz) do j=1+jlow, ny-jup flux(:,1) = (fxi(1:nvars,nx+1,j,nz,1)-fxi(1:nvars,nx,j,nz,1))/dx flux(:,2) = (fxi(1:nvars,nx,j+1,nz,2)-fxi(1:nvars,nx,j,nz,2))/dy flux(:,3) = (fxi(1:nvars,nx,j,nz+1,3)-fxi(1:nvars,nx,j,nz,3))/dz uxloc = ux(:,nx,j,nz) vxloc = vx(:,nx,j,nz) uxref = (ux(:,nx,j,nz+1)+ux(:,nx+1,j,nz))*0.5d0 p = (vx(5,nx,j,nz+1)+vx(5,nx+1,j,nz))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,j,nz), & flux, Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,j,nz), sat, & flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, j, nz) ! rotation back call InverseRotateXYZ(source(:,nx,j,nz), sat, flux, Proj, & cbc_ixup) ! Correct flux to the last point source(:,nx,j,nz) = source(:,nx,j,nz) - sat(:) & - flux(:,1) - flux(:,2) -flux(:,3) enddo endif !----------------------------------------------------------- ! ! BACK BOTTOM EDGE ! !----------------------------------------------------------- if ( cbc_iylow .eq. cbc_izlow .and. cbc_iylow .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_iylow ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_iylow ) iup = 1 ! low priority corner (just avoid) if ( cbc_iylow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 endif Proj = 0.0d0 Proj(1,2) = 1.0d0 Proj(2,1) = -isqr2 Proj(2,3) = -isqr2 Proj(3,1) = -isqr2 Proj(3,3) = isqr2 dl = sqrt(dy*dz) do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,1,1,1)-fxi(1:nvars,i,1,1,1))/dx flux(:,2) = (fxi(1:nvars,i,2,1,2) -fxi(1:nvars,i,1,1,2))/dy flux(:,3) = (fxi(1:nvars,i,1,2,3) -fxi(1:nvars,i,1,1,3))/dz uxloc = ux(:,i,1,1) vxloc = vx(:,i,1,1) uxref = (ux(:,i,1,0)+ux(:,i,0,1))*0.5d0 p = (vx(5,i,1,0)+vx(5,i,0,1))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,1,1), flux, & Proj, cbc_iylow) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,1,1), sat, & flux, p, dl, cbc_force_matrix(:,3), cbc_iylow, i, 1, 1) ! rotation back call InverseRotateXYZ(source(:,i,1,1), sat, flux, Proj, cbc_iylow) ! Correct flux to the last point source(:,i,1,1) = source(:,i,1,1) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) enddo endif !----------------------------------------------------------- ! ! FRONT BOTTOM EDGE ! !----------------------------------------------------------- if ( cbc_iylow .eq. cbc_izup .and. cbc_iylow .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_iylow ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_iylow ) iup = 1 ! low priority corner (just avoid) if ( cbc_iylow .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 endif Proj = 0.0d0 Proj(1,2) = 1.0d0 Proj(2,1) = -isqr2 Proj(2,3) = isqr2 Proj(3,1) = isqr2 Proj(3,3) = isqr2 dl = sqrt(dy*dz) do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,1,nz,1)-fxi(1:nvars,i,1,nz,1))/dx flux(:,2) = (fxi(1:nvars,i,2,nz,2) -fxi(1:nvars,i,1,nz,2))/dy flux(:,3) = (fxi(1:nvars,i,1,nz+1,3)-fxi(1:nvars,i,1,nz,3))/dz uxloc = ux(:,i,1,nz) vxloc = vx(:,i,1,nz) uxref = (ux(:,i,1,nz+1)+ux(:,i,0,nz))*0.5d0 p = (vx(5,i,1,nz+1)+vx(5,i,0,nz))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,1,nz), flux,& Proj, cbc_iylow) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,1,nz), sat, & flux, p, dl, cbc_force_matrix(:,3), cbc_iylow, i, 1, nz) ! rotation back call InverseRotateXYZ(source(:,i,1,nz), sat, flux, Proj, & cbc_iylow) ! Correct flux to the last point source(:,i,1,nz) = source(:,i,1,nz) - sat(:) & - flux(:,1) -flux(:,2) -flux(:,3) enddo endif !----------------------------------------------------------- ! ! BACK TOP EDGE ! !----------------------------------------------------------- if ( cbc_iyup .eq. cbc_izlow .and. cbc_iyup .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_iyup ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_iyup ) iup = 1 ! low priority corner (just avoid) if ( cbc_iyup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 endif Proj = 0.0d0 Proj(1,2) = -1.0d0 Proj(2,1) = isqr2 Proj(2,3) = isqr2 Proj(3,1) = -isqr2 Proj(3,3) = isqr2 dl = sqrt(dy*dz) do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,ny,1,1)-fxi(1:nvars,i,ny,1,1))/dx flux(:,2) = (fxi(1:nvars,i,ny+1,1,2)-fxi(1:nvars,i,ny,1,2))/dy flux(:,3) = (fxi(1:nvars,i,ny,2,3) -fxi(1:nvars,i,ny,1,3))/dz uxloc = ux(:,i,ny,1) vxloc = vx(:,i,ny,1) uxref = (ux(:,i,ny,0)+ux(:,i,ny+1,1))*0.5d0 p = (vx(5,i,ny,0)+vx(5,i,ny+1,1))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,ny,1), flux,& Proj, cbc_iyup) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,ny,1), sat, & flux, p, dl, cbc_force_matrix(:,4), cbc_iyup, i, ny, 1) ! rotation back call InverseRotateXYZ(source(:,i,ny,1), sat, flux, Proj, cbc_iyup) ! Correct flux to the last point source(:,i,ny,1) = source(:,i,ny,1) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) enddo endif !----------------------------------------------------------- ! ! FRONT TOP EDGE ! !----------------------------------------------------------- if ( cbc_iyup .eq. cbc_izup .and. cbc_iyup .ne. CLES_CBC_NONE ) then ilow = 0; iup = 0 ! lower corner if ( cbc_ixlow .eq. cbc_iyup ) ilow = 1 ! upper corner if ( cbc_ixup .eq. cbc_iyup ) iup = 1 ! low priority corner (just avoid) if ( cbc_iyup .eq. CLES_CBC_OUTFLOW ) then if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1 if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1 endif Proj = 0.0d0 Proj(1,2) = -1.0d0 Proj(2,1) = isqr2 Proj(2,3) = -isqr2 Proj(3,1) = isqr2 Proj(3,3) = isqr2 dl = sqrt(dy*dz) do i=1+ilow, nx-iup flux(:,1) = (fxi(1:nvars,i+1,ny,nz,1)-fxi(1:nvars,i,ny,nz,1))/dx flux(:,2) = (fxi(1:nvars,i,ny+1,nz,2)-fxi(1:nvars,i,ny,nz,2))/dy flux(:,3) = (fxi(1:nvars,i,ny,nz+1,3)-fxi(1:nvars,i,ny,nz,3))/dz uxloc = ux(:,i,ny,nz) vxloc = vx(:,i,ny,nz) uxref = (ux(:,i,ny,nz+1)+ux(:,i,ny+1,nz))*0.5d0 p = (vx(5,i,ny,nz+1)+vx(5,i,ny+1,nz))*0.5d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,ny,nz), & flux, Proj, cbc_iyup) call SingleCSAT(uxloc, vxloc, uxref, source(:,i,ny,nz), sat, & flux, p, dl, cbc_force_matrix(:,4), cbc_iyup, i, ny, nz) ! rotation back call InverseRotateXYZ(source(:,i,ny,nz), sat, flux, Proj, & cbc_iyup) ! Correct flux to the last point source(:,i,ny,nz) = source(:,i,ny,nz) - sat(:) & - flux(:,1) - flux(:,2) -flux(:,3) enddo endif !----------------------------------------------------------- ! ! LEFT-BOTTOM-BACK CORNER ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_iylow .and. cbc_ixlow .ne. CLES_CBC_NONE & .and. cbc_iylow .eq. cbc_izlow ) then call setproj(-3.0d0*pi/4.0d0,-phi, Proj) flux(:,1) = (fxi(1:nvars,2,1,1,1)-fxi(1:nvars,1,1,1,1))/dx flux(:,2) = (fxi(1:nvars,1,2,1,2)-fxi(1:nvars,1,1,1,2))/dy flux(:,3) = (fxi(1:nvars,1,1,2,3)-fxi(1:nvars,1,1,1,3))/dz uxloc = ux(:,1,1,1) vxloc = vx(:,1,1,1) uxref = (ux(:,0,1,1)+ux(:,1,0,1)+ux(:,1,1,0))/3.0d0 dl = (dx*dy*dz)**(1.0d0/3.0d0) p = (vx(5,0,1,1)+vx(5,1,0,1)+vx(5,1,1,0))/3.0d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,1,1), flux, & Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,1,1), sat, & flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, 1, 1) ! rotation back call InverseRotateXYZ(source(:,1,1,1), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,1,1) = source(:,1,1,1) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) endif !----------------------------------------------------------- ! ! RIGHT-BOTTOM-BACK CORNER ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_iylow .and. cbc_ixup .ne. CLES_CBC_NONE & .and. cbc_iylow .eq. cbc_izlow ) then call setproj(-pi/4.0d0,-phi, Proj) flux(:,1) = (fxi(1:nvars,nx+1,1,1,1)-fxi(1:nvars,nx,1,1,1))/dx flux(:,2) = (fxi(1:nvars,nx,2,1,2) -fxi(1:nvars,nx,1,1,2))/dy flux(:,3) = (fxi(1:nvars,nx,1,2,3) -fxi(1:nvars,nx,1,1,3))/dz uxloc = ux(:,nx,1,1) vxloc = vx(:,nx,1,1) uxref = (ux(:,nx+1,1,1)+ux(:,nx,0,1)+ux(:,nx,1,0))/3.0d0 dl = (dx*dy*dz)**(1.0d0/3.0d0) p = (vx(5,nx+1,1,1)+vx(5,nx,0,1)+vx(5,nx,1,0))/3.0d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,1,1), flux, & Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,1,1), sat, & flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, 1, 1) ! rotation back call InverseRotateXYZ(source(:,nx,1,1), sat, flux, Proj, cbc_ixup) ! Correct flux to the last point source(:,nx,1,1) = source(:,nx,1,1) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) endif !----------------------------------------------------------- ! ! LEFT-TOP-BACK CORNER ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_iyup .and. cbc_ixlow .ne. CLES_CBC_NONE & .and. cbc_iyup .eq. cbc_izlow ) then call setproj(3.0d0*pi/4.0d0,-phi, Proj) flux(:,1) = (fxi(1:nvars,2,ny,1,1) -fxi(1:nvars,1,ny,1,1))/dx flux(:,2) = (fxi(1:nvars,1,ny+1,1,2)-fxi(1:nvars,1,ny,1,2))/dy flux(:,3) = (fxi(1:nvars,1,ny,2,3) -fxi(1:nvars,1,ny,1,3))/dz uxloc = ux(:,1,ny,1) vxloc = vx(:,1,ny,1) uxref = (ux(:,0,ny,1)+ux(:,1,ny+1,1)+ux(:,1,ny,0))/3.0d0 dl = (dx*dy*dz)**(1.0d0/3.0d0) p = (vx(5,0,ny,1)+vx(5,1,ny+1,1)+vx(5,1,ny,0))/3.0d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,ny,1), flux, & Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,ny,1), sat, & flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, ny, 1) ! rotation back call InverseRotateXYZ(source(:,1,ny,1), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,ny,1) = source(:,1,ny,1) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) endif !----------------------------------------------------------- ! ! LEFT-BOTTOM-FRONT CORNER ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_iylow .and. cbc_ixlow .ne. CLES_CBC_NONE & .and. cbc_iylow .eq. cbc_izup ) then call setproj(-3.0d0*pi/4.0d0,phi, Proj) flux(:,1) = (fxi(1:nvars,2,1,nz,1) -fxi(1:nvars,1,1,nz,1))/dx flux(:,2) = (fxi(1:nvars,1,2,nz,2) -fxi(1:nvars,1,1,nz,2))/dy flux(:,3) = (fxi(1:nvars,1,1,nz+1,3)-fxi(1:nvars,1,1,nz,3))/dz uxloc = ux(:,1,1,nz) vxloc = vx(:,1,1,nz) uxref = (ux(:,0,1,nz)+ux(:,1,0,nz)+ux(:,1,1,nz+1))/3.0d0 dl = (dx*dy*dz)**(1.0d0/3.0d0) p = (vx(5,0,1,nz)+vx(5,1,0,nz)+vx(5,1,1,nz+1))/3.0d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,1,nz), flux, & Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,1,nz), sat, & flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, 1, nz) ! rotation back call InverseRotateXYZ(source(:,1,1,nz), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,1,nz) = source(:,1,1,nz) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) endif !----------------------------------------------------------- ! ! RIGHT-TOP-BACK CORNER ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_iyup .and. cbc_ixup .ne. CLES_CBC_NONE & .and. cbc_iyup .eq. cbc_izlow ) then call setproj(pi/4.0d0,-phi, Proj) flux(:,1) = (fxi(1:nvars,nx+1,ny,1,1)-fxi(1:nvars,nx,ny,1,1))/dx flux(:,2) = (fxi(1:nvars,nx,ny+1,1,2)-fxi(1:nvars,nx,ny,1,2))/dy flux(:,3) = (fxi(1:nvars,nx,ny,2,3) -fxi(1:nvars,nx,ny,1,3))/dz uxloc = ux(:,nx,ny,1) vxloc = vx(:,nx,ny,1) uxref = (ux(:,nx+1,ny,1)+ux(:,nx,ny+1,1)+ux(:,nx,ny,0))/3.0d0 dl = (dx*dy*dz)**(1.0d0/3.0d0) p = (vx(5,nx+1,ny,1)+vx(5,nx,ny+1,1)+vx(5,nx,ny,0))/3.0d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,ny,1), flux, & Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,ny,1), sat, & flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, ny, 1) ! rotation back call InverseRotateXYZ(source(:,nx,ny,1), sat, flux, Proj, cbc_ixup) ! Correct flux to the last point source(:,nx,ny,1) = source(:,nx,ny,1) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) endif !----------------------------------------------------------- ! ! LEFT-TOP-FRONT CORNER ! !----------------------------------------------------------- if ( cbc_ixlow .eq. cbc_iyup .and. cbc_ixlow .ne. CLES_CBC_NONE & .and. cbc_iyup .eq. cbc_izup ) then call setproj(3.0d0*pi/4.0d0,phi, Proj) flux(:,1) = (fxi(1:nvars,2,ny,nz,1) -fxi(1:nvars,1,ny,nz,1))/dx flux(:,2) = (fxi(1:nvars,1,ny+1,nz,2)-fxi(1:nvars,1,ny,nz,2))/dy flux(:,3) = (fxi(1:nvars,1,ny,nz+1,3)-fxi(1:nvars,1,ny,nz,3))/dz uxloc = ux(:,1,ny,nz) vxloc = vx(:,1,ny,nz) uxref = (ux(:,0,ny,nz)+ux(:,1,ny+1,nz)+ux(:,1,ny,nz+1))/3.0d0 dl = (dx*dy*dz)**(1.0d0/3.0d0) p = (vx(5,0,ny,nz)+vx(5,1,ny+1,nz)+vx(5,1,ny,nz+1))/3.0d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,ny,nz), flux, & Proj, cbc_ixlow) call SingleCSAT(uxloc, vxloc, uxref, source(:,1,ny,nz), sat, & flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, ny, nz) ! rotation back call InverseRotateXYZ(source(:,1,ny,nz), sat, flux, Proj, cbc_ixlow) ! Correct flux to the last point source(:,1,ny,nz) = source(:,1,ny,nz) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) endif !----------------------------------------------------------- ! ! RIGHT-BOTTOM-FRONT CORNER ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_iylow .and. cbc_ixup .ne. CLES_CBC_NONE & .and. cbc_iylow .eq. cbc_izup ) then call setproj(-pi/4.0d0,phi, Proj) flux(:,1) = (fxi(1:nvars,nx+1,1,nz,1)-fxi(1:nvars,nx,1,nz,1))/dx flux(:,2) = (fxi(1:nvars,nx,2,nz,2) -fxi(1:nvars,nx,1,nz,2))/dy flux(:,3) = (fxi(1:nvars,nx,1,nz+1,3)-fxi(1:nvars,nx,1,nz,3))/dz uxloc = ux(:,nx,1,nz) vxloc = vx(:,nx,1,nz) uxref = (ux(:,nx+1,1,nz)+ux(:,nx,0,nz)+ux(:,nx,1,nz+1))/3.0d0 dl = (dx*dy*dz)**(1.0d0/3.0d0) p = (vx(5,nx+1,1,nz)+vx(5,nx,0,nz)+vx(5,nx,1,nz+1))/3.0d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,1,nz), flux, & Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,1,nz), sat, & flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, 1, nz) ! rotation back call InverseRotateXYZ(source(:,nx,1,nz), sat, flux, Proj, cbc_ixup) ! Correct flux to the last point source(:,nx,1,nz) = source(:,nx,1,nz) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) endif !----------------------------------------------------------- ! ! RIGHT-TOP-FRONT CORNER ! !----------------------------------------------------------- if ( cbc_ixup .eq. cbc_iyup .and. cbc_ixup .ne. CLES_CBC_NONE & .and. cbc_iyup .eq. cbc_izup ) then call setproj(pi/4.0d0,phi, Proj) flux(:,1) = (fxi(1:nvars,nx+1,ny,nz,1)-fxi(1:nvars,nx,ny,nz,1))/dx flux(:,2) = (fxi(1:nvars,nx,ny+1,nz,2)-fxi(1:nvars,nx,ny,nz,2))/dy flux(:,3) = (fxi(1:nvars,nx,ny,nz+1,3)-fxi(1:nvars,nx,ny,nz,3))/dz uxloc = ux(:,nx,ny,nz) vxloc = vx(:,nx,ny,nz) uxref = (ux(:,nx+1,ny,nz)+ux(:,nx,ny+1,nz)+ux(:,nx,ny,nz+1))/3.0d0 dl = (dx*dy*dz)**(1.0d0/3.0d0) p = (vx(5,nx+1,ny,nz)+vx(5,nx,ny+1,nz)+vx(5,nx,ny,nz+1))/3.0d0 ! rotate call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,ny,nz), flux, & Proj, cbc_ixup) call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,ny,nz), sat, & flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, ny, nz) ! rotation back call InverseRotateXYZ(source(:,nx,ny,nz), sat, flux, Proj, cbc_ixup) ! Correct flux to the last point source(:,nx,ny,nz) = source(:,nx,ny,nz) - sat(:) & - flux(:,1) -flux(:,2) - flux(:,3) endif RETURN END SUBROUTINE ThreeDAcousticBC SUBROUTINE OneDAcousticViscousBC(is, fivx, direction) ! ---- Shared Variables USE mesh USE array_bounds USE method_parms USE Interp_coeffs USE Generic_EvalGamma IMPLICIT NONE include 'cles.i' DOUBLE PRECISION, INTENT(INOUT) :: fivx(ixlo:ixhi) INTEGER, INTENT(IN) :: direction, is logical :: c1, c2 IF ( cbc_active .NE. CLES_TRUE ) RETURN IF ( direction .EQ. CLES_CBC_XDIR ) THEN IF ( cbc_ixlow .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on left-most cell wall ! so that in this boundary cell ! df/dx = ( fivx(1) - fivx(2) )/dx =0 c1 = (cbc_ixlow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_ixlow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(1) = fivx(2) ENDIF IF ( cbc_ixup .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on right-most cell wall ! so that in this boundary cell ! df/dx = ( fivx(nx+1) - fivx(nx) )/dx =0 c1 = (cbc_ixup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_ixup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(nx+1) = fivx(nx) ENDIF ENDIF RETURN END SUBROUTINE OneDAcousticViscousBC SUBROUTINE TwoDAcousticViscousBC(is, fivx, direction) ! ---- Shared Variables USE mesh USE array_bounds USE method_parms USE Interp_coeffs USE Generic_EvalGamma IMPLICIT NONE include 'cles.i' DOUBLE PRECISION, INTENT(INOUT) :: fivx(ixlo:ixhi,iylo:iyhi) INTEGER, INTENT(IN) :: direction, is logical :: c1, c2 IF ( cbc_active .NE. CLES_TRUE ) RETURN IF ( direction .EQ. CLES_CBC_XDIR ) THEN IF ( cbc_ixlow .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on left-most cell wall ! so that in this boundary cell ! df/dx = ( fivx(1,:) - fivx(2,:) )/dx =0 c1 = (cbc_ixlow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_ixlow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(1,:) = fivx(2,:) ENDIF IF ( cbc_ixup .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on right-most cell wall ! so that in this boundary cell ! df/dx = ( fivx(nx+1,:) - fivx(nx,:) )/dx =0 c1 = (cbc_ixup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_ixup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(nx+1,:) = fivx(nx,:) ENDIF ENDIF IF ( direction .EQ. CLES_CBC_YDIR ) THEN IF ( cbc_iylow .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on left-most cell wall ! so that in this boundary cell ! df/dy = ( fivx(:,1) - fivx(:,2) )/dy =0 c1 = (cbc_iylow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_iylow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(:,1) = fivx(:,2) ENDIF IF ( cbc_iyup .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on right-most cell wall ! so that in this boundary cell ! df/dy = ( fivx(:,ny+1) - fivx(:,ny) )/dy =0 c1 = (cbc_iyup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_iyup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(:,ny+1) = fivx(:,ny) ENDIF ENDIF RETURN END SUBROUTINE TwoDAcousticViscousBC SUBROUTINE ThreeDAcousticViscousBC(is, fivx, direction) ! ---- Shared Variables USE mesh USE array_bounds USE method_parms USE Interp_coeffs USE Generic_EvalGamma IMPLICIT NONE include 'cles.i' DOUBLE PRECISION, INTENT(INOUT) :: fivx(ixlo:ixhi,iylo:iyhi,izlo:izhi) INTEGER, INTENT(IN) :: direction, is logical :: c1, c2 IF ( cbc_active .NE. CLES_TRUE ) RETURN IF ( direction .EQ. CLES_CBC_XDIR ) THEN IF ( cbc_ixlow .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on left-most cell wall ! so that in this boundary cell ! df/dx = ( fivx(1,:,:) - fivx(2,:,:) )/dx =0 c1 = (cbc_ixlow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_ixlow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(1,:,:) = fivx(2,:,:) ENDIF IF ( cbc_ixup .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on right-most cell wall ! so that in this boundary cell ! df/dx = ( fivx(nx+1,:,:) - fivx(nx,:,:) )/dx =0 c1 = (cbc_ixup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_ixup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(nx+1,:,:) = fivx(nx,:,:) ENDIF ENDIF IF ( direction .EQ. CLES_CBC_YDIR ) THEN IF ( cbc_iylow .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on left-most cell wall ! so that in this boundary cell ! df/dy = ( fivx(:,1,:) - fivx(:,2,:) )/dy =0 c1 = (cbc_iylow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_iylow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(:,1,:) = fivx(:,2,:) ENDIF IF ( cbc_iyup .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on right-most cell wall ! so that in this boundary cell ! df/dy = ( fivx(:,ny+1,:) - fivx(:,ny,:) )/dy =0 c1 = (cbc_iyup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_iyup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(:,ny+1,:) = fivx(:,ny,:) ENDIF ENDIF IF ( direction .EQ. CLES_CBC_ZDIR ) THEN IF ( cbc_izlow .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on left-most cell wall ! so that in this boundary cell ! df/dz = ( fivx(:,:,1) - fivx(:,:,2) )/dz =0 c1 = (cbc_izlow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_izlow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(:,:,1) = fivx(:,:,2) ENDIF IF ( cbc_izup .NE. CLES_CBC_NONE ) THEN ! Inflow ! set the flux on right-most cell wall ! so that in this boundary cell ! df/dz = ( fivx(:,:,nz+1) - fivx(:,:,nz) )/dz =0 c1 = (cbc_izup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction) c2 = (cbc_izup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction) if ( c1 .or. c2 ) fivx(:,:,nz+1) = fivx(:,:,nz) ENDIF ENDIF RETURN END SUBROUTINE ThreeDAcousticViscousBC END MODULE Generic_AcousticBC SUBROUTINE DontUseAcousticBC(ival) use mesh USE Generic_AcousticBC USE method_parms IMPLICIT NONE include 'cles.i' INTEGER, INTENT(IN) :: ival integer iret cbc_active = CLES_FALSE if ( ALLOCATED(cbc_force_matrix) ) then DEALLOCATE(cbc_force_matrix, stat=iret) if ( iret .ne. 0 ) then call cleslog_error(CLESLOG_ERROR_DEALLOCATE, & 'DontUseAcousticBC/cbc_force_matrix', iret) endif endif RETURN END SUBROUTINE DontUseAcousticBC SUBROUTINE UseAcousticBC(ival) use mesh USE Generic_AcousticBC USE method_parms IMPLICIT NONE include 'cles.i' INTEGER, INTENT(IN) :: ival integer iret cbc_active = CLES_TRUE cbc_mode = ival ALLOCATE(cbc_force_matrix(nvars,6), stat=iret) if ( iret .ne. 0 ) then call cleslog_error(CLESLOG_ERROR_ALLOCATE, & 'UseAcousticBC/cbc_force_matrix', iret) endif cbc_force_matrix(:,:) = 0.0d0 RETURN END SUBROUTINE UseAcousticBC SUBROUTINE SetAcousticBoundary(idir, iside, itype) use method_parms IMPLICIT NONE include 'cles.i' INTEGER, INTENT(IN) :: idir, iside, itype if ( idir .lt. 1 .or. idir .gt. 3 ) then print *, 'AcousticBC Error: boundary direction incorrect' stop endif if ( iside .lt. CLES_CBC_LEFT .or. iside .gt. CLES_CBC_RIGHT ) then print *, 'AcousticBC Error: boundary side incorrect' stop endif cbc_direction(iside + (idir-1)*2) = itype END SUBROUTINE SetAcousticBoundary SUBROUTINE SetAcousticForceMatrix(idir, iside, vec) USE Generic_AcousticBC USE method_parms IMPLICIT NONE include 'cles.i' INTEGER, INTENT(IN) :: idir, iside DOUBLE PRECISION, INTENT(IN) :: vec(nvars) if ( idir .lt. 1 .or. idir .gt. 3 ) then print *, 'AcousticBC Error: boundary direction incorrect' stop endif if ( iside .lt. CLES_CBC_LEFT .or. iside .gt. CLES_CBC_RIGHT ) then print *, 'AcousticBC Error: boundary side incorrect' stop endif cbc_force_matrix(:,iside + (idir-1)*2) = vec(:) END SUBROUTINE SetAcousticForceMatrix FUNCTION CBCIsActive() USE Generic_AcousticBC IMPLICIT NONE include 'cles.i' LOGICAL :: CBCIsActive if ( cbc_active .eq. CLES_TRUE ) then CBCIsActive = .TRUE. else CBCIsActive = .FALSE. endif RETURN END FUNCTION CBCIsActive SUBROUTINE SetAcousticImpedance(sigma) USE Generic_AcousticBC IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: sigma cbc_sigma = sigma cbc_factor = cbc_sigma*(1.0d0-cbc_Mach*cbc_Mach)*0.5d0 END SUBROUTINE SetAcousticImpedance SUBROUTINE SetAcousticMach(Ma1) USE Generic_AcousticBC IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: Ma1 cbc_Mach = Ma1 cbc_factor = cbc_sigma*(1.0d0-cbc_Mach*cbc_Mach)*0.5d0 END SUBROUTINE SetAcousticMach