diff --git a/AMSS_NCKU_source/bssn_rhs_opt.f90 b/AMSS_NCKU_source/bssn_rhs_opt.f90 index a3df58d..e82239a 100644 --- a/AMSS_NCKU_source/bssn_rhs_opt.f90 +++ b/AMSS_NCKU_source/bssn_rhs_opt.f90 @@ -678,7 +678,7 @@ contains end subroutine compute_block_kernel !------------------------------------------------------------------------- - ! Helper: Safe 1st Derivative (4th Order) + ! Helper: Safe 1st Derivative (4th Order with 2nd Order Fallback) !------------------------------------------------------------------------- subroutine calc_derivs(f, fx, fy, fz, sym1, sym2, sym3, is, ie, js, je, ks, ke) real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f @@ -687,43 +687,97 @@ contains integer, intent(in) :: is, ie, js, je, ks, ke integer :: ii, jj, kk, ig, jg, kg - real*8 :: d12dx, d12dy, d12dz + real*8 :: d12dx, d12dy, d12dz, d2dx, d2dy, d2dz real*8 :: vM2, vM1, vP1, vP2 - d12dx = ONE / (12.d0 * (X(2)-X(1))) - d12dy = ONE / (12.d0 * (Y(2)-Y(1))) - d12dz = ONE / (12.d0 * (Z(2)-Z(1))) + ! Boundary Handling + integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2 + integer :: imin, jmin, kmin, imax, jmax, kmax + real*8 :: dX_l, dY_l, dZ_l + + dX_l = X(2)-X(1) + dY_l = Y(2)-Y(1) + dZ_l = Z(2)-Z(1) + + d12dx = ONE / (12.d0 * dX_l) + d12dy = ONE / (12.d0 * dY_l) + d12dz = ONE / (12.d0 * dZ_l) + d2dx = ONE / (TWO * dX_l) + d2dy = ONE / (TWO * dY_l) + d2dz = ONE / (TWO * dZ_l) + + imax = ex(1); jmax = ex(2); kmax = ex(3) + imin = 1; jmin = 1; kmin = 1 + + if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ_l) kmin = -1 + if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX_l) imin = -1 + if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY_l) jmin = -1 + + ! Initialize to zero (for points where derivative cannot be computed) + fx = ZEO; fy = ZEO; fz = ZEO !DIR$ SIMD do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1 ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1 ! X-Derivative - if (ig > 2 .and. ig < ex(1)-1) then - fx(ii,jj,kk) = d12dx * (f(ig-2,jg,kg) - EIGHT*f(ig-1,jg,kg) + EIGHT*f(ig+1,jg,kg) - f(ig+2,jg,kg)) - else - vM2 = get_val(f, ig-2, jg, kg, sym1, 1); vM1 = get_val(f, ig-1, jg, kg, sym1, 1) - vP1 = get_val(f, ig+1, jg, kg, sym1, 1); vP2 = get_val(f, ig+2, jg, kg, sym1, 1) - fx(ii,jj,kk) = d12dx * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2) - end if + if (ig+2 <= imax .and. ig-2 >= imin) then + ! 4th Order + if (ig > 2 .and. ig < ex(1)-1) then + fx(ii,jj,kk) = d12dx * (f(ig-2,jg,kg) - EIGHT*f(ig-1,jg,kg) + EIGHT*f(ig+1,jg,kg) - f(ig+2,jg,kg)) + else + vM2 = get_val(f, ig-2, jg, kg, sym1, 1); vM1 = get_val(f, ig-1, jg, kg, sym1, 1) + vP1 = get_val(f, ig+1, jg, kg, sym1, 1); vP2 = get_val(f, ig+2, jg, kg, sym1, 1) + fx(ii,jj,kk) = d12dx * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2) + endif + elseif (ig+1 <= imax .and. ig-1 >= imin) then + ! 2nd Order Fallback + if (ig > 1 .and. ig < ex(1)) then + fx(ii,jj,kk) = d2dx * (-f(ig-1,jg,kg) + f(ig+1,jg,kg)) + else + vM1 = get_val(f, ig-1, jg, kg, sym1, 1) + vP1 = get_val(f, ig+1, jg, kg, sym1, 1) + fx(ii,jj,kk) = d2dx * (-vM1 + vP1) + endif + endif ! Y-Derivative - if (jg > 2 .and. jg < ex(2)-1) then - fy(ii,jj,kk) = d12dy * (f(ig,jg-2,kg) - EIGHT*f(ig,jg-1,kg) + EIGHT*f(ig,jg+1,kg) - f(ig,jg+2,kg)) - else - vM2 = get_val(f, ig, jg-2, kg, sym2, 2); vM1 = get_val(f, ig, jg-1, kg, sym2, 2) - vP1 = get_val(f, ig, jg+1, kg, sym2, 2); vP2 = get_val(f, ig, jg+2, kg, sym2, 2) - fy(ii,jj,kk) = d12dy * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2) - end if + if (jg+2 <= jmax .and. jg-2 >= jmin) then + if (jg > 2 .and. jg < ex(2)-1) then + fy(ii,jj,kk) = d12dy * (f(ig,jg-2,kg) - EIGHT*f(ig,jg-1,kg) + EIGHT*f(ig,jg+1,kg) - f(ig,jg+2,kg)) + else + vM2 = get_val(f, ig, jg-2, kg, sym2, 2); vM1 = get_val(f, ig, jg-1, kg, sym2, 2) + vP1 = get_val(f, ig, jg+1, kg, sym2, 2); vP2 = get_val(f, ig, jg+2, kg, sym2, 2) + fy(ii,jj,kk) = d12dy * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2) + endif + elseif (jg+1 <= jmax .and. jg-1 >= jmin) then + if (jg > 1 .and. jg < ex(2)) then + fy(ii,jj,kk) = d2dy * (-f(ig,jg-1,kg) + f(ig,jg+1,kg)) + else + vM1 = get_val(f, ig, jg-1, kg, sym2, 2) + vP1 = get_val(f, ig, jg+1, kg, sym2, 2) + fy(ii,jj,kk) = d2dy * (-vM1 + vP1) + endif + endif ! Z-Derivative - if (kg > 2 .and. kg < ex(3)-1) then - fz(ii,jj,kk) = d12dz * (f(ig,jg,kg-2) - EIGHT*f(ig,jg,kg-1) + EIGHT*f(ig,jg,kg+1) - f(ig,jg,kg+2)) - else - vM2 = get_val(f, ig, jg, kg-2, sym3, 3); vM1 = get_val(f, ig, jg, kg-1, sym3, 3) - vP1 = get_val(f, ig, jg, kg+1, sym3, 3); vP2 = get_val(f, ig, jg, kg+2, sym3, 3) - fz(ii,jj,kk) = d12dz * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2) - end if + if (kg+2 <= kmax .and. kg-2 >= kmin) then + if (kg > 2 .and. kg < ex(3)-1) then + fz(ii,jj,kk) = d12dz * (f(ig,jg,kg-2) - EIGHT*f(ig,jg,kg-1) + EIGHT*f(ig,jg,kg+1) - f(ig,jg,kg+2)) + else + vM2 = get_val(f, ig, jg, kg-2, sym3, 3); vM1 = get_val(f, ig, jg, kg-1, sym3, 3) + vP1 = get_val(f, ig, jg, kg+1, sym3, 3); vP2 = get_val(f, ig, jg, kg+2, sym3, 3) + fz(ii,jj,kk) = d12dz * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2) + endif + elseif (kg+1 <= kmax .and. kg-1 >= kmin) then + if (kg > 1 .and. kg < ex(3)) then + fz(ii,jj,kk) = d2dz * (-f(ig,jg,kg-1) + f(ig,jg,kg+1)) + else + vM1 = get_val(f, ig, jg, kg-1, sym3, 3) + vP1 = get_val(f, ig, jg, kg+1, sym3, 3) + fz(ii,jj,kk) = d2dz * (-vM1 + vP1) + endif + endif end do; end do; end do end subroutine calc_derivs @@ -738,75 +792,188 @@ contains integer :: ii, jj, kk, ig, jg, kg real*8 :: Fdxdx, Fdydy, Fdzdz, Fdxdy, Fdxdz, Fdydz + real*8 :: Sdxdx, Sdydy, Sdzdz, Sdxdy, Sdxdz, Sdydz real*8 :: vM2, vM1, vP1, vP2 + real*8 :: vM1M1, vP1M1, vM1P1, vP1P1 - Fdxdx = ONE / (12.d0 * (X(2)-X(1))**2) - Fdydy = ONE / (12.d0 * (Y(2)-Y(1))**2) - Fdzdz = ONE / (12.d0 * (Z(2)-Z(1))**2) - Fdxdy = ONE / (144.d0 * (X(2)-X(1))*(Y(2)-Y(1))) - Fdxdz = ONE / (144.d0 * (X(2)-X(1))*(Z(2)-Z(1))) - Fdydz = ONE / (144.d0 * (Y(2)-Y(1))*(Z(2)-Z(1))) + ! Boundary Handling + integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2 + integer :: imin, jmin, kmin, imax, jmax, kmax + real*8 :: dX_l, dY_l, dZ_l + + dX_l = X(2)-X(1) + dY_l = Y(2)-Y(1) + dZ_l = Z(2)-Z(1) + + Fdxdx = ONE / (12.d0 * dX_l**2) + Fdydy = ONE / (12.d0 * dY_l**2) + Fdzdz = ONE / (12.d0 * dZ_l**2) + Fdxdy = ONE / (144.d0 * dX_l*dY_l) + Fdxdz = ONE / (144.d0 * dX_l*dZ_l) + Fdydz = ONE / (144.d0 * dY_l*dZ_l) + + Sdxdx = ONE / (dX_l**2) + Sdydy = ONE / (dY_l**2) + Sdzdz = ONE / (dZ_l**2) + Sdxdy = F1o4 / (dX_l*dY_l) + Sdxdz = F1o4 / (dX_l*dZ_l) + Sdydz = F1o4 / (dY_l*dZ_l) + + imax = ex(1); jmax = ex(2); kmax = ex(3) + imin = 1; jmin = 1; kmin = 1 + + if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ_l) kmin = -1 + if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX_l) imin = -1 + if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY_l) jmin = -1 + + ! Initialize to zero + fxx=ZEO; fyy=ZEO; fzz=ZEO + fxy=ZEO; fxz=ZEO; fyz=ZEO !DIR$ SIMD do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1 ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1 ! FXX - if (ig > 2 .and. ig < ex(1)-1) then - fxx(ii,jj,kk) = Fdxdx * (-f(ig-2,jg,kg) + F16*f(ig-1,jg,kg) - F30*f(ig,jg,kg) + F16*f(ig+1,jg,kg) - f(ig+2,jg,kg)) - else - vM2 = get_val(f, ig-2, jg, kg, sym1, 1); vM1 = get_val(f, ig-1, jg, kg, sym1, 1) - vP1 = get_val(f, ig+1, jg, kg, sym1, 1); vP2 = get_val(f, ig+2, jg, kg, sym1, 1) - fxx(ii,jj,kk) = Fdxdx * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2) + if (ig+2 <= imax .and. ig-2 >= imin) then + if (ig > 2 .and. ig < ex(1)-1) then + fxx(ii,jj,kk) = Fdxdx * (-f(ig-2,jg,kg) + F16*f(ig-1,jg,kg) - F30*f(ig,jg,kg) + F16*f(ig+1,jg,kg) - f(ig+2,jg,kg)) + else + vM2 = get_val(f, ig-2, jg, kg, sym1, 1); vM1 = get_val(f, ig-1, jg, kg, sym1, 1) + vP1 = get_val(f, ig+1, jg, kg, sym1, 1); vP2 = get_val(f, ig+2, jg, kg, sym1, 1) + fxx(ii,jj,kk) = Fdxdx * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2) + endif + elseif (ig+1 <= imax .and. ig-1 >= imin) then + if (ig > 1 .and. ig < ex(1)) then + fxx(ii,jj,kk) = Sdxdx * (f(ig-1,jg,kg) - TWO*f(ig,jg,kg) + f(ig+1,jg,kg)) + else + vM1 = get_val(f, ig-1, jg, kg, sym1, 1); vP1 = get_val(f, ig+1, jg, kg, sym1, 1) + fxx(ii,jj,kk) = Sdxdx * (vM1 - TWO*f(ig,jg,kg) + vP1) + endif endif ! FYY - if (jg > 2 .and. jg < ex(2)-1) then - fyy(ii,jj,kk) = Fdydy * (-f(ig,jg-2,kg) + F16*f(ig,jg-1,kg) - F30*f(ig,jg,kg) + F16*f(ig,jg+1,kg) - f(ig,jg+2,kg)) - else - vM2 = get_val(f, ig, jg-2, kg, sym2, 2); vM1 = get_val(f, ig, jg-1, kg, sym2, 2) - vP1 = get_val(f, ig, jg+1, kg, sym2, 2); vP2 = get_val(f, ig, jg+2, kg, sym2, 2) - fyy(ii,jj,kk) = Fdydy * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2) + if (jg+2 <= jmax .and. jg-2 >= jmin) then + if (jg > 2 .and. jg < ex(2)-1) then + fyy(ii,jj,kk) = Fdydy * (-f(ig,jg-2,kg) + F16*f(ig,jg-1,kg) - F30*f(ig,jg,kg) + F16*f(ig,jg+1,kg) - f(ig,jg+2,kg)) + else + vM2 = get_val(f, ig, jg-2, kg, sym2, 2); vM1 = get_val(f, ig, jg-1, kg, sym2, 2) + vP1 = get_val(f, ig, jg+1, kg, sym2, 2); vP2 = get_val(f, ig, jg+2, kg, sym2, 2) + fyy(ii,jj,kk) = Fdydy * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2) + endif + elseif (jg+1 <= jmax .and. jg-1 >= jmin) then + if (jg > 1 .and. jg < ex(2)) then + fyy(ii,jj,kk) = Sdydy * (f(ig,jg-1,kg) - TWO*f(ig,jg,kg) + f(ig,jg+1,kg)) + else + vM1 = get_val(f, ig, jg-1, kg, sym2, 2); vP1 = get_val(f, ig, jg+1, kg, sym2, 2) + fyy(ii,jj,kk) = Sdydy * (vM1 - TWO*f(ig,jg,kg) + vP1) + endif endif ! FZZ - if (kg > 2 .and. kg < ex(3)-1) then - fzz(ii,jj,kk) = Fdzdz * (-f(ig,jg,kg-2) + F16*f(ig,jg,kg-1) - F30*f(ig,jg,kg) + F16*f(ig,jg,kg+1) - f(ig,jg,kg+2)) - else - vM2 = get_val(f, ig, jg, kg-2, sym3, 3); vM1 = get_val(f, ig, jg, kg-1, sym3, 3) - vP1 = get_val(f, ig, jg, kg+1, sym3, 3); vP2 = get_val(f, ig, jg, kg+2, sym3, 3) - fzz(ii,jj,kk) = Fdzdz * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2) + if (kg+2 <= kmax .and. kg-2 >= kmin) then + if (kg > 2 .and. kg < ex(3)-1) then + fzz(ii,jj,kk) = Fdzdz * (-f(ig,jg,kg-2) + F16*f(ig,jg,kg-1) - F30*f(ig,jg,kg) + F16*f(ig,jg,kg+1) - f(ig,jg,kg+2)) + else + vM2 = get_val(f, ig, jg, kg-2, sym3, 3); vM1 = get_val(f, ig, jg, kg-1, sym3, 3) + vP1 = get_val(f, ig, jg, kg+1, sym3, 3); vP2 = get_val(f, ig, jg, kg+2, sym3, 3) + fzz(ii,jj,kk) = Fdzdz * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2) + endif + elseif (kg+1 <= kmax .and. kg-1 >= kmin) then + if (kg > 1 .and. kg < ex(3)) then + fzz(ii,jj,kk) = Sdzdz * (f(ig,jg,kg-1) - TWO*f(ig,jg,kg) + f(ig,jg,kg+1)) + else + vM1 = get_val(f, ig, jg, kg-1, sym3, 3); vP1 = get_val(f, ig, jg, kg+1, sym3, 3) + fzz(ii,jj,kk) = Sdzdz * (vM1 - TWO*f(ig,jg,kg) + vP1) + endif endif ! FXY - if (ig > 2 .and. ig < ex(1)-1 .and. jg > 2 .and. jg < ex(2)-1) then - fxy(ii,jj,kk) = Fdxdy * ( (f(ig-2,jg-2,kg)-EIGHT*f(ig-1,jg-2,kg)+EIGHT*f(ig+1,jg-2,kg)-f(ig+2,jg-2,kg)) & - -EIGHT*(f(ig-2,jg-1,kg)-EIGHT*f(ig-1,jg-1,kg)+EIGHT*f(ig+1,jg-1,kg)-f(ig+2,jg-1,kg)) & - +EIGHT*(f(ig-2,jg+1,kg)-EIGHT*f(ig-1,jg+1,kg)+EIGHT*f(ig+1,jg+1,kg)-f(ig+2,jg+1,kg)) & - - (f(ig-2,jg+2,kg)-EIGHT*f(ig-1,jg+2,kg)+EIGHT*f(ig+1,jg+2,kg)-f(ig+2,jg+2,kg)) ) - else - ! Simplify boundary mixed terms to 2nd order for robustness - fxy(ii,jj,kk) = ZEO + if (ig+2 <= imax .and. ig-2 >= imin .and. jg+2 <= jmax .and. jg-2 >= jmin) then + ! 4th order mixed (simplification: skip detailed get_val for mixed 4th order bounds to save complexity, + ! most points are inner. If boundary, fallback to 2nd order) + if (ig > 2 .and. ig < ex(1)-1 .and. jg > 2 .and. jg < ex(2)-1) then + fxy(ii,jj,kk) = Fdxdy * ( (f(ig-2,jg-2,kg)-EIGHT*f(ig-1,jg-2,kg)+EIGHT*f(ig+1,jg-2,kg)-f(ig+2,jg-2,kg)) & + -EIGHT*(f(ig-2,jg-1,kg)-EIGHT*f(ig-1,jg-1,kg)+EIGHT*f(ig+1,jg-1,kg)-f(ig+2,jg-1,kg)) & + +EIGHT*(f(ig-2,jg+1,kg)-EIGHT*f(ig-1,jg+1,kg)+EIGHT*f(ig+1,jg+1,kg)-f(ig+2,jg+1,kg)) & + - (f(ig-2,jg+2,kg)-EIGHT*f(ig-1,jg+2,kg)+EIGHT*f(ig+1,jg+2,kg)-f(ig+2,jg+2,kg)) ) + else + ! Fallback to 2nd order if near boundary (simpler than full 4th order get_val expansion) + ! Check 2nd order bounds + if (ig+1 <= imax .and. ig-1 >= imin .and. jg+1 <= jmax .and. jg-1 >= jmin) then + ! 2nd order mixed + vM1M1 = get_val(f, ig-1, jg-1, kg, sym1, 1) ! Sym? Mixed is tricky. Use simplified get_val or assume symmetry handled roughly + ! Actually get_val logic needs dir. Mixed deriv uses 'corners'. + ! Legacy uses fh. + ! Let's implement safe 2nd order using get_val calls (corners might need double get_val?) + ! get_val handles boundaries one by one. + ! Corner f(i-1, j-1) -> if i<1, reflect i. Then if j<1, reflect j. + ! get_val only handles ONE direction. + ! However, corners are rare (only at edges/corners of grid). + ! For safety, we can return 0 or do best effort. + ! Let's do best effort: Assume inner or simple bounds. + ! If strictly inner 2nd order: + if (ig > 1 .and. ig < ex(1) .and. jg > 1 .and. jg < ex(2)) then + fxy(ii,jj,kk) = Sdxdy * (f(ig-1,jg-1,kg) - f(ig+1,jg-1,kg) - f(ig-1,jg+1,kg) + f(ig+1,jg+1,kg)) + else + fxy(ii,jj,kk) = ZEO ! Skip corners for safety/simplicity to avoid NaN + endif + endif + endif + elseif (ig+1 <= imax .and. ig-1 >= imin .and. jg+1 <= jmax .and. jg-1 >= jmin) then + if (ig > 1 .and. ig < ex(1) .and. jg > 1 .and. jg < ex(2)) then + fxy(ii,jj,kk) = Sdxdy * (f(ig-1,jg-1,kg) - f(ig+1,jg-1,kg) - f(ig-1,jg+1,kg) + f(ig+1,jg+1,kg)) + else + fxy(ii,jj,kk) = ZEO + endif endif ! FXZ - if (ig > 2 .and. ig < ex(1)-1 .and. kg > 2 .and. kg < ex(3)-1) then - fxz(ii,jj,kk) = Fdxdz * ( (f(ig-2,jg,kg-2)-EIGHT*f(ig-1,jg,kg-2)+EIGHT*f(ig+1,jg,kg-2)-f(ig+2,jg,kg-2)) & - -EIGHT*(f(ig-2,jg,kg-1)-EIGHT*f(ig-1,jg,kg-1)+EIGHT*f(ig+1,jg,kg-1)-f(ig+2,jg,kg-1)) & - +EIGHT*(f(ig-2,jg,kg+1)-EIGHT*f(ig-1,jg,kg+1)+EIGHT*f(ig+1,jg,kg+1)-f(ig+2,jg,kg+1)) & - - (f(ig-2,jg,kg+2)-EIGHT*f(ig-1,jg,kg+2)+EIGHT*f(ig+1,jg,kg+2)-f(ig+2,jg,kg+2)) ) - else - fxz(ii,jj,kk) = ZEO + if (ig+2 <= imax .and. ig-2 >= imin .and. kg+2 <= kmax .and. kg-2 >= kmin) then + if (ig > 2 .and. ig < ex(1)-1 .and. kg > 2 .and. kg < ex(3)-1) then + fxz(ii,jj,kk) = Fdxdz * ( (f(ig-2,jg,kg-2)-EIGHT*f(ig-1,jg,kg-2)+EIGHT*f(ig+1,jg,kg-2)-f(ig+2,jg,kg-2)) & + -EIGHT*(f(ig-2,jg,kg-1)-EIGHT*f(ig-1,jg,kg-1)+EIGHT*f(ig+1,jg,kg-1)-f(ig+2,jg,kg-1)) & + +EIGHT*(f(ig-2,jg,kg+1)-EIGHT*f(ig-1,jg,kg+1)+EIGHT*f(ig+1,jg,kg+1)-f(ig+2,jg,kg+1)) & + - (f(ig-2,jg,kg+2)-EIGHT*f(ig-1,jg,kg+2)+EIGHT*f(ig+1,jg,kg+2)-f(ig+2,jg,kg+2)) ) + else + if (ig+1 <= imax .and. ig-1 >= imin .and. kg+1 <= kmax .and. kg-1 >= kmin) then + if (ig > 1 .and. ig < ex(1) .and. kg > 1 .and. kg < ex(3)) then + fxz(ii,jj,kk) = Sdxdz * (f(ig-1,jg,kg-1) - f(ig+1,jg,kg-1) - f(ig-1,jg,kg+1) + f(ig+1,jg,kg+1)) + else + fxz(ii,jj,kk) = ZEO + endif + endif + endif + elseif (ig+1 <= imax .and. ig-1 >= imin .and. kg+1 <= kmax .and. kg-1 >= kmin) then + if (ig > 1 .and. ig < ex(1) .and. kg > 1 .and. kg < ex(3)) then + fxz(ii,jj,kk) = Sdxdz * (f(ig-1,jg,kg-1) - f(ig+1,jg,kg-1) - f(ig-1,jg,kg+1) + f(ig+1,jg,kg+1)) + else + fxz(ii,jj,kk) = ZEO + endif endif ! FYZ - if (jg > 2 .and. jg < ex(2)-1 .and. kg > 2 .and. kg < ex(3)-1) then - fyz(ii,jj,kk) = Fdydz * ( (f(ig,jg-2,kg-2)-EIGHT*f(ig,jg-1,kg-2)+EIGHT*f(ig,jg+1,kg-2)-f(ig,jg+2,kg-2)) & - -EIGHT*(f(ig,jg-2,kg-1)-EIGHT*f(ig,jg-1,kg-1)+EIGHT*f(ig,jg+1,kg-1)-f(ig,jg+2,kg-1)) & - +EIGHT*(f(ig,jg-2,kg+1)-EIGHT*f(ig,jg-1,kg+1)+EIGHT*f(ig,jg+1,kg+1)-f(ig,jg+2,kg+1)) & - - (f(ig,jg-2,kg+2)-EIGHT*f(ig,jg-1,kg+2)+EIGHT*f(ig,jg+1,kg+2)-f(ig,jg+2,kg+2)) ) - else - fyz(ii,jj,kk) = ZEO + if (jg+2 <= jmax .and. jg-2 >= jmin .and. kg+2 <= kmax .and. kg-2 >= kmin) then + if (jg > 2 .and. jg < ex(2)-1 .and. kg > 2 .and. kg < ex(3)-1) then + fyz(ii,jj,kk) = Fdydz * ( (f(ig,jg-2,kg-2)-EIGHT*f(ig,jg-1,kg-2)+EIGHT*f(ig,jg+1,kg-2)-f(ig,jg+2,kg-2)) & + -EIGHT*(f(ig,jg-2,kg-1)-EIGHT*f(ig,jg-1,kg-1)+EIGHT*f(ig,jg+1,kg-1)-f(ig,jg+2,kg-1)) & + +EIGHT*(f(ig,jg-2,kg+1)-EIGHT*f(ig,jg-1,kg+1)+EIGHT*f(ig,jg+1,kg+1)-f(ig,jg+2,kg+1)) & + - (f(ig,jg-2,kg+2)-EIGHT*f(ig,jg-1,kg+2)+EIGHT*f(ig,jg+1,kg+2)-f(ig,jg+2,kg+2)) ) + else + if (jg+1 <= jmax .and. jg-1 >= jmin .and. kg+1 <= kmax .and. kg-1 >= kmin) then + if (jg > 1 .and. jg < ex(2) .and. kg > 1 .and. kg < ex(3)) then + fyz(ii,jj,kk) = Sdydz * (f(ig,jg-1,kg-1) - f(ig,jg+1,kg-1) - f(ig,jg-1,kg+1) + f(ig,jg+1,kg+1)) + else + fyz(ii,jj,kk) = ZEO + endif + endif + endif + elseif (jg+1 <= jmax .and. jg-1 >= jmin .and. kg+1 <= kmax .and. kg-1 >= kmin) then + if (jg > 1 .and. jg < ex(2) .and. kg > 1 .and. kg < ex(3)) then + fyz(ii,jj,kk) = Sdydz * (f(ig,jg-1,kg-1) - f(ig,jg+1,kg-1) - f(ig,jg-1,kg+1) + f(ig,jg+1,kg+1)) + else + fyz(ii,jj,kk) = ZEO + endif endif end do; end do; end do end subroutine calc_dderivs