Homec4science
Diffusion miniapp_gbs_gradients (mathieu_bench_omp_cuda_m100) Browse

R11986/grad_fd4.h
f16fd9b8f371mathieu_bench_omp_cuda_…

Last Change
Raw File

grad_fd4.h
View Options

!****************************************************
!************** FIRST DERIVATIVE ********************
!****************************************************
!-----------------------------------------------------------------------------
subroutine gradz_n2n_fd4(f , f_z)
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_z
integer :: iz
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:) = deltazi*coef_der1_n2n(:)
f_z(:,:,:) = nan_
do iz=izs,ize
f_z(:, :, iz) = coef_der(1)*f(:, :, iz-2) &
+ coef_der(2)*f(:, :, iz-1) &
+ coef_der(3)*f(:, :, iz+1) &
+ coef_der(4)*f(:, :, iz+2)
end do
end subroutine gradz_n2n_fd4
!----------------------------------------------------------------------------
subroutine grady_n2n_fd4(f,f_y)
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_y
integer :: iy
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:)= deltayi*coef_der1_n2n(:)
f_y (:,:,:) = nan_
do iy=iys,iye
f_y(iy, :, :) = coef_der(1)*f(iy-2, :, :) &
+ coef_der(2)*f(iy-1, :, :) &
+ coef_der(3)*f(iy+1, :, :) &
+ coef_der(4)*f(iy+2, :, :)
end do
end subroutine grady_n2n_fd4
!---------------------------------------------------------------------------
subroutine gradx_n2n_fd4(f,f_x)
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_x
integer :: ix
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:)= deltaxi*coef_der1_n2n(:)
f_x (:,:,: ) = nan_
do ix=ixs,ixe
f_x(:, ix, :) = coef_der(1)*f(:, ix-2, :) &
+ coef_der(2)*f(:, ix-1, :) &
+ coef_der(3)*f(:, ix+1, :) &
+ coef_der(4)*f(:, ix+2, :)
end do
end subroutine gradx_n2n_fd4
!----------------------------------------------------------------------------
subroutine gradz_v2n_fd4(f,f_z_v2n)
use prec_const
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_z_v2n
integer :: iy, iz
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:) = deltazi*coef_der1_stag(:)
f_z_v2n(:,:,:) = nan_
do iy=iys,iye
do iz=izs,ize
f_z_v2n(iy,:,iz) = coef_int(1)*( coef_der(1)*f(iy-1, :, iz-1) &
+ coef_der(2)*f(iy-1, :, iz) &
+ coef_der(3)*f(iy-1, :, iz +1) &
+ coef_der(4)*f(iy-1, :, iz+2) )&
+ coef_int(2)*( coef_der(1)*f(iy, :, iz-1) &
+ coef_der(2)*f(iy, :, iz) &
+ coef_der(3)*f(iy, :, iz+1) &
+ coef_der(4)*f(iy, :, iz+2) )&
+ coef_int(3)*( coef_der(1)*f(iy+1, :, iz-1) &
+ coef_der(2)*f(iy+1, :, iz) &
+ coef_der(3)*f(iy+1, :, iz+1) &
+ coef_der(4)*f(iy+1, :, iz+2) )&
+ coef_int(4)*( coef_der(1)*f(iy+2, :, iz-1) &
+ coef_der(2)*f(iy+2, :, iz) &
+ coef_der(3)*f(iy+2, :, iz+1) &
+ coef_der(4)*f(iy+2, :, iz+2) )
end do
end do
end subroutine gradz_v2n_fd4
!-------------------------------------------------------------------------------------------
subroutine grady_v2n_fd4(f,f_y_v2n)
use prec_const
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_y_v2n
integer :: iy, iz
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:) = deltayi*coef_der1_stag(:)
f_y_v2n(:,:,:) = nan_
do iy=iys,iye
do iz=izs,ize
f_y_v2n(iy,:,iz) = coef_int(1)*( coef_der(1)*f(iy-1, :, iz-1) &
+ coef_der(2)*f(iy, :, iz-1) &
+ coef_der(3)*f(iy+1, :, iz-1) &
+ coef_der(4)*f(iy+2, :, iz-1) )&
+ coef_int(2)*( coef_der(1)*f(iy-1, :, iz) &
+ coef_der(2)*f(iy, :, iz) &
+ coef_der(3)*f(iy+1, :, iz) &
+ coef_der(4)*f(iy+2, :, iz) )&
+ coef_int(3)*( coef_der(1)*f(iy-1, :, iz+1) &
+ coef_der(2)*f(iy, :, iz+1) &
+ coef_der(3)*f(iy+1, :, iz+1) &
+ coef_der(4)*f(iy+2, :, iz+1) )&
+ coef_int(4)*( coef_der(1)*f(iy-1, :, iz+2) &
+ coef_der(2)*f(iy, :, iz+2) &
+ coef_der(3)*f(iy+1, :, iz+2) &
+ coef_der(4)*f(iy+2, :, iz+2) )
end do
end do
end subroutine grady_v2n_fd4
!-----------------------------------------------------------------------------------------
subroutine gradz_n2v_fd4(f,f_z_n2v)
use prec_const
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_z_n2v
integer :: iy, iz
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:) = deltazi*coef_der1_stag(:)
f_z_n2v(:,:,:) = nan_
do iy=iys,iye
do iz=izs,ize
f_z_n2v(iy,:,iz) = coef_int(1)*( coef_der(1)*f(iy-2, :, iz-2) &
+ coef_der(2)*f(iy-2, :, iz-1) &
+ coef_der(3)*f(iy-2, :, iz ) &
+ coef_der(4)*f(iy-2, :, iz+1) )&
+ coef_int(2)*( coef_der(1)*f(iy-1, :, iz-2) &
+ coef_der(2)*f(iy-1, :, iz-1) &
+ coef_der(3)*f(iy-1, :, iz) &
+ coef_der(4)*f(iy-1, :, iz+1) )&
+ coef_int(3)*( coef_der(1)*f(iy, :, iz-2) &
+ coef_der(2)*f(iy, :, iz-1) &
+ coef_der(3)*f(iy, :, iz) &
+ coef_der(4)*f(iy, :, iz+1) )&
+ coef_int(4)*( coef_der(1)*f(iy+1, :, iz-2) &
+ coef_der(2)*f(iy+1, :, iz-1) &
+ coef_der(3)*f(iy+1, :, iz) &
+ coef_der(4)*f(iy+1, :, iz+1) )
end do
end do
end subroutine gradz_n2v_fd4
!-----------------------------------------------------------------------------------
subroutine grady_n2v_fd4(f,f_y_n2v)
use prec_const
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_y_n2v
integer :: iy, iz
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:) = deltayi*coef_der1_stag(:)
f_y_n2v(:,:,:) = nan_
do iy=iys,iye
do iz=izs,ize
f_y_n2v(iy,:,iz) = coef_int(1)*( coef_der(1)*f(iy-2, :, iz-2) &
+ coef_der(2)*f(iy-1, :, iz-2) &
+ coef_der(3)*f(iy, :, iz-2) &
+ coef_der(4)*f(iy+1, :, iz-2) )&
+ coef_int(2)*( coef_der(1)*f(iy-2, :, iz-1) &
+ coef_der(2)*f(iy-1, :, iz-1) &
+ coef_der(3)*f(iy, :, iz-1) &
+ coef_der(4)*f(iy+1, :, iz-1) )&
+ coef_int(3)*( coef_der(1)*f(iy-2, :, iz) &
+ coef_der(2)*f(iy-1, :, iz) &
+ coef_der(3)*f(iy, :, iz) &
+ coef_der(4)*f(iy+1, :, iz) )&
+ coef_int(4)*( coef_der(1)*f(iy-2, :, iz+1) &
+ coef_der(2)*f(iy-1, :, iz+1) &
+ coef_der(3)*f(iy, :, iz+1) &
+ coef_der(4)*f(iy+1, :, iz+1) )
end do
end do
end subroutine grady_n2v_fd4
!-----------------------------------------------------------------------------------
!***************************************************
!*************** SECOND DERIVATIVE *****************
!***************************************************
subroutine gradxx_n2n_fd4( f , f_xx )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg),intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_xx
real(dp), DIMENSION(1:5) :: coef_der
integer ::ix
coef_der(:) = (deltaxisq)*coef_der2(:)
f_xx(:,:,:) = nan_
do ix=ixs,ixe
f_xx(:, ix, :) = coef_der(1)*f(:, ix-2, :) &
+ coef_der(2)*f(:, ix-1, :) &
+ coef_der(3)*f(:, ix , :) &
+ coef_der(4)*f(:, ix+1, :) &
+ coef_der(5)*f(:, ix+2, :)
end do
end subroutine gradxx_n2n_fd4
!-----------------------------------------------------------------------------------
subroutine gradzz_n2n_fd4(f,f_zz)
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_zz
integer :: iz
real(dp), DIMENSION(1:5) :: coef_der
coef_der(:) = (deltazisq)*coef_der2(:)
f_zz(:, :,:) = nan_
do iz=izs,ize
f_zz(:, :, iz) = coef_der(1)*f(:, :, iz-2) &
+ coef_der(2)*f(:, :, iz-1) &
+ coef_der(3)*f(:, :, iz) &
+ coef_der(4)*f(:, :, iz+1) &
+ coef_der(5)*f(:, :, iz+2)
end do
end subroutine gradzz_n2n_fd4
!-------------------------------------------------------------------------------------
subroutine gradyy_n2n_fd4(f,f_yy)
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_yy
integer :: iy
real(dp), DIMENSION(1:5) :: coef_der
coef_der(:) = deltayisq*coef_der2(:)
f_yy(:,:,:) = nan_
do iy=iys,iye
f_yy(iy, :, :) = coef_der(1)*f(iy-2, :, :) &
+ coef_der(2)*f(iy-1, :, :) &
+ coef_der(3)*f(iy, :, :) &
+ coef_der(4)*f(iy+1, :, :) &
+ coef_der(5)*f(iy+2, :, :)
end do
end subroutine gradyy_n2n_fd4
subroutine gradyy_v2n_fd4(f,f_yy_v2n)
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg),intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_yy_v2n
integer :: iy,iz
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:) = deltayisq*coef_der2_stag(:)
!this derivative is only at 2nd order. to do it at 4th we would have needed 3 ghost cells
f_yy_v2n (:,:,:) = nan_
do iy= iys, iye
do iz = izs,ize
f_yy_v2n (iy,:,iz) = coef_int(1)*( coef_der(1)*f(iy-1,:,iz-1) &
+ coef_der(2)*f(iy ,:,iz-1) &
+ coef_der(3)*f(iy+1,:,iz-1) &
+ coef_der(4)*f(iy+2,:,iz-1))&
+ coef_int(2)* ( coef_der(1)*f(iy-1,:,iz) &
+ coef_der(2)*f(iy ,:,iz) &
+ coef_der(3)*f(iy+1,:,iz) &
+ coef_der(4)*f(iy+2,:,iz))&
+ coef_int(3)*( coef_der(1)*f(iy-1,:,iz+1) &
+ coef_der(2)*f(iy ,:,iz+1) &
+ coef_der(3)*f(iy+1,:,iz+1) &
+ coef_der(4)*f(iy+2,:,iz+1))&
+ coef_int(4)*( coef_der(1)*f(iy-1,:,iz+2) &
+ coef_der(2)*f(iy ,:,iz+2) &
+ coef_der(3)*f(iy+1,:,iz+2) &
+ coef_der(4)*f(iy+2,:,iz+2))
end do
end do
end subroutine gradyy_v2n_fd4
subroutine gradyy_n2v_fd4 (f, f_yy_n2v)
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg),intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_yy_n2v
integer :: iy,iz
real(dp), DIMENSION(1:4) :: coef_der
coef_der(:) = deltayisq*coef_der2_stag(:)
f_yy_n2v (:,:,:) = nan_
!this derivative is only at 2nd order. to do it at 4th we would have needed 3 ghost cells
do iy= iys, iye
do iz = izs, ize
f_yy_n2v (iy,:,iz) = coef_int(1)*( coef_der(1)*f(iy-2 ,:,iz-2) &
+ coef_der(2)*f(iy-1 ,:,iz-2) &
+ coef_der(3)*f(iy ,:,iz-2) &
+ coef_der(4)*f(iy+1,:,iz-2) ) &
+ coef_int(2)*( coef_der(1)*f(iy-2 ,:,iz-1) &
+ coef_der(2)*f(iy-1 ,:,iz-1) &
+ coef_der(3)*f(iy ,:,iz-1) &
+ coef_der(4)*f(iy+1 ,:,iz-1) ) &
+ coef_int(3)*( coef_der(1)*f(iy-2 ,:,iz ) &
+ coef_der(2)*f(iy-1,:,iz ) &
+ coef_der(3)*f(iy ,:,iz ) &
+ coef_der(4)*f(iy+1,:,iz ) ) &
+ coef_int(4)*( coef_der(1)*f(iy-2 ,:,iz+1) &
+ coef_der(2)*f(iy-1 ,:,iz+1) &
+ coef_der(3)*f(iy ,:,iz+1) &
+ coef_der(4)*f(iy+1 ,:,iz+1) )
end do
end do
end subroutine gradyy_n2v_fd4
!-------------------------------------------------------------------------------------
!******************************************************
!***************** MIXED DERIVATIVES ******************
!******************************************************
subroutine gradyz_n2v_fd4(f, f_yz_n2v)
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg),intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_yz_n2v
integer :: iy, iz
real(dp), DIMENSION(1:4) :: coef_der_y, coef_der_z
coef_der_z(:) = deltazi*coef_der1_stag(:)
coef_der_y(:) = deltayi*coef_der1_stag(:)
f_yz_n2v(:,:,:) = nan_
do iy=iys,iye
do iz=izs,ize
f_yz_n2v(iy,:,iz) = coef_der_y(1)*( coef_der_z(1)*f(iy-2, :, iz-2) &
+ coef_der_z(2)*f(iy-2, :, iz-1) &
+ coef_der_z(3)*f(iy-2, :, iz ) &
+ coef_der_z(4)*f(iy-2, :, iz+1) )&
+ coef_der_y(2)*( coef_der_z(1)*f(iy-1, :, iz-2) &
+ coef_der_z(2)*f(iy-1, :, iz-1) &
+ coef_der_z(3)*f(iy-1, :, iz) &
+ coef_der_z(4)*f(iy-1, :, iz+1) )&
+ coef_der_y(3)*( coef_der_z(1)*f(iy, :, iz-2) &
+ coef_der_z(2)*f(iy, :, iz-1) &
+ coef_der_z(3)*f(iy, :, iz) &
+ coef_der_z(4)*f(iy, :, iz+1) )&
+ coef_der_y(4)*( coef_der_z(1)*f(iy+1, :, iz-2) &
+ coef_der_z(2)*f(iy+1, :, iz-1) &
+ coef_der_z(3)*f(iy+1, :, iz) &
+ coef_der_z(4)*f(iy+1, :, iz+1) )
end do
end do
end subroutine gradyz_n2v_fd4
!-------------------------------------------------------------------------------------
subroutine gradyz_v2n_fd4(f, f_yz_v2n)
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg),intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_yz_v2n
real(dp), DIMENSION(1:4) :: coef_der_y, coef_der_z
integer :: iy, iz
coef_der_z(:) = deltazi*coef_der1_stag(:)
coef_der_y(:) = deltayi*coef_der1_stag(:)
f_yz_v2n(:,:,:) = nan_
do iy=iys,iye
do iz=izs,ize
f_yz_v2n(iy,:,iz) = coef_der_y(1)*( coef_der_z(1)*f(iy-1, :, iz-1) &
+ coef_der_z(2)*f(iy-1, :, iz) &
+ coef_der_z(3)*f(iy-1, :, iz +1) &
+ coef_der_z(4)*f(iy-1, :, iz+2) )&
+ coef_der_y(2)*( coef_der_z(1)*f(iy, :, iz-1) &
+ coef_der_z(2)*f(iy, :, iz) &
+ coef_der_z(3)*f(iy, :, iz+1) &
+ coef_der_z(4)*f(iy, :, iz+2) )&
+ coef_der_y(3)*( coef_der_z(1)*f(iy+1, :, iz-1) &
+ coef_der_z(2)*f(iy+1, :, iz) &
+ coef_der_z(3)*f(iy+1, :, iz+1) &
+ coef_der_z(4)*f(iy+1, :, iz+2) )&
+ coef_der_y(4)*( coef_der_z(1)*f(iy+2, :, iz-1) &
+ coef_der_z(2)*f(iy+2, :, iz) &
+ coef_der_z(3)*f(iy+2, :, iz+1) &
+ coef_der_z(4)*f(iy+2, :, iz+2) )
end do
end do
end subroutine gradyz_v2n_fd4
!---------------------------------------------------------------------------------
subroutine scalyy_fd4( f, diff, fyy )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f, diff
real(dp), dimension(iysg:iyeg ,ixsg:ixeg ,izsg:izeg), intent(out) :: fyy
real(dp), dimension(iysg:iyeg ,ixsg:ixeg ,izsg:izeg) :: grady_diff, grady_f
integer :: iy
real(dp), DIMENSION(1:5) :: coef_der
coef_der(:) = deltayisq*coef_der2(:)
fyy(:, :, :) = nan_
call grady_n2n_fd4(f, grady_f)
call grady_n2n_fd4(diff, grady_diff)
do iy=iys,iye
fyy(iy, :, :) = diff(iy,:,:)*( coef_der(1)*f(iy-2, :, :) &
+ coef_der(2)*f(iy-1, :, :) &
+ coef_der(3)*f(iy , :, :) &
+ coef_der(4)*f(iy+1, :, :) &
+ coef_der(5)*f(iy+2, :, :) ) &
+ grady_diff(iy,:,:)*grady_f(iy,:,:)
end do
end subroutine scalyy_fd4
!--------------------------------------------------------------------------------
!Diffusion operators
subroutine scalxx_fd4( f, diff, fxx )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f, diff
real(dp), dimension(iysg:iyeg ,ixsg:ixeg ,izsg:izeg), intent(out) :: fxx
real(dp), dimension(iysg:iyeg ,ixsg:ixeg ,izsg:izeg) :: gradx_diff, gradx_f
integer :: ix
real(dp), DIMENSION(1:5) :: coef_der
coef_der(:) = deltaxisq*coef_der2(:)
fxx(:, :, :) = nan_
call gradx_n2n_fd4(f, gradx_f)
call gradx_n2n_fd4(diff, gradx_diff)
do ix=ixs,ixe
fxx(:, ix, :) = diff(:,ix,:)*( coef_der(1)*f(:, ix-2, :) &
+ coef_der(2)*f(:, ix-1, :) &
+ coef_der(3)*f(:, ix , :) &
+ coef_der(4)*f(:, ix+1, :) &
+ coef_der(5)*f(:, ix+2, :) ) &
+ gradx_diff(:,ix,:)*gradx_f(:,ix,:)
end do
end subroutine scalxx_fd4
!--------------------------------------------------------------------
subroutine lapl_fd4(f, f_lapl)
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg),intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: f_lapl
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_xx, f_yy
integer:: ix,iy,iz
f_lapl(:,:,:) = nan_
call gradxx_n2n_fd4(f, f_xx)
call gradyy_n2n_fd4(f, f_yy)
do iz = izs,ize
do ix = ixs,ixe
do iy = iys,iye
f_lapl(iy,ix,iz) = f_xx(iy,ix,iz) + f_yy(iy,ix,iz)
end do
end do
end do
end subroutine lapl_fd4
!*********************************************************
!******************* PARALLEL GRADIENT *******************
!*********************************************************
! parallel gradient for finite differences 4rth order from v grid to n grid
subroutine gradpar_v2n_fd4 ( f, flu, frd , f_grad)
!flu and frd are useless here for periodic BC
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_grad
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_z, f_y, f_n, f_x_n !f_x_v
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2), intent(in) :: flu, frd
integer :: ix, iy
f_grad(:,:,:) = nan_
call gradz_v2n_fd4(f, f_z)
call grady_v2n_fd4(f, f_y)
! its important to do the int first
call interp_v2n_fd4(f, flu, frd, f_n )
call gradx_n2n_fd4(f_n, f_x_n)
!call gradx_n2n_fd4(f, f_x_n)
do ix = ixs,ixe
do iy = iys, iye
f_grad(iy,ix,izs:ize) = gradpar_z*f_z(iy,ix,izs:ize) + gradpar_y_n(iy,ix)*f_y(iy,ix,izs:ize) &
+ gradpar_x_n(iy,ix)*f_x_n(iy,ix,izs:ize)
end do
end do
end subroutine gradpar_v2n_fd4
!----------------------------------------------------------------------------
! parallel gradient for finite differences 4rth order from n grid to v grid
subroutine gradpar_n2v_fd4(f, flu, frd , f_grad)
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) , intent(out) :: f_grad
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_z,f_y,f_v,f_x_v !f_x_n
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2), intent(in) :: flu, frd
integer:: ix, iy
f_grad(:,:,:)=nan_
call gradz_n2v_fd4(f, f_z)
call grady_n2v_fd4(f,f_y)
! its important to do the int first
call interp_n2v_fd4(f, flu, frd,f_v)
call gradx_n2n_fd4(f_v, f_x_v)
do ix= ixs,ixe
do iy= iys, iye
!f_grad(iy,ix,:) = gradpar_z*f_z(iy,ix,:) + gradpar_y_v(iy,ix)*f_y(iy,ix,:) &
! + gradpar_x_v(iy,ix)*f_x_v(iy,ix,:)
f_grad(iy,ix,:) = f_v(iy,ix,:)
end do
end do
end subroutine gradpar_n2v_fd4
!----------------------------------------------------------------------------
! parallel gradient for finite differences 4rth order from n grid to n grid
subroutine gradpar_n2n_fd4(f, flu, frd , f_grad)
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: f_grad
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_z,f_y,f_x
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2), intent(in) :: flu, frd
integer :: ix, iy
f_grad(:,:,:) = nan_
call gradz_n2n_fd4(f, f_z)
call grady_n2n_fd4(f, f_y)
call gradx_n2n_fd4(f, f_x)
do ix = ixs, ixe
do iy = iys, iye
f_grad(iy,ix,izs:ize) = gradpar_z*f_z(iy,ix,izs:ize) + gradpar_y_n(iy,ix)*f_y(iy,ix,izs:ize) &
+ gradpar_x_n(iy,ix)*f_x(iy,ix,izs:ize)
end do
end do
end subroutine gradpar_n2n_fd4
!-------------------------------------------------------------------------------
!----------------------------------------------------------------------------
! parallel gradient for finite differences 4rth order from n grid to n grid
subroutine gradpar_v2v_fd4(f, flu, frd, f_grad)
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: f_grad
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_z,f_y,f_x
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2), intent(in) :: flu, frd
integer :: ix, iy
f_grad(:,:,:) = nan_
call gradz_n2n_fd4(f,f_z )
call grady_n2n_fd4(f,f_y)
call gradx_n2n_fd4(f,f_x)
do ix = ixs, ixe
do iy = iys, iye
f_grad(iy,ix,izs:ize) = gradpar_z*f_z(iy,ix,izs:ize) + gradpar_y_v(iy,ix)*f_y(iy,ix,izs:ize)&
+ gradpar_x_v(iy,ix)*f_x(iy,ix,izs:ize)
end do
end do
end subroutine gradpar_v2v_fd4
!-------------------------------------------------------------------------------
subroutine grad2par_n2n_fd4(f,flu,frd, f_lapl)
! flu and frd are useless here bum passing them anyway to respect the code structure in numerics
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2), intent(in) :: flu, frd
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_lapl
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_zz, f_yy, f_xx
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_zy, f_yz, f_zx, f_xz, f_xy, f_yx
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_x, f_y, f_z
integer :: ix,iy
call gradzz_n2n_fd4(f ,f_zz)
call gradyy_n2n_fd4(f ,f_yy)
call gradxx_n2n_fd4(f ,f_xx)
call gradz_n2n_fd4 (f ,f_z )
call gradx_n2n_fd4 (f ,f_x )
call grady_n2n_fd4 (f ,f_y )
call grady_n2n_fd4 (f_x,f_yx)
call gradx_n2n_fd4 (f_y,f_xy)
call gradx_n2n_fd4 (f_z,f_xz)
call gradz_n2n_fd4 (f_x,f_zx)
call grady_n2n_fd4 (f_z,f_yz)
call gradz_n2n_fd4 (f_y,f_zy)
f_lapl(:,:,:) = nan_
do ix = ixs,ixe
do iy = iys, iye
f_lapl(iy,ix,izs:ize) = f_zz(iy,ix,izs:ize) &
+ grad2par_yy_n(iy,ix) *f_yy(iy,ix,izs:ize) &
+ grad2par_xx_n(iy,ix) *f_xx(iy,ix,izs:ize) &
+ grad2par_zy_n(iy,ix) *(f_zy(iy,ix,izs:ize)+f_yz(iy,ix,izs:ize)) &
+ grad2par_zx_n(iy,ix) *(f_zx(iy,ix,izs:ize)+f_xz(iy,ix,izs:ize)) &
+ grad2par_xy_n(iy,ix) *(f_xy(iy,ix,izs:ize)+f_yx(iy,ix,izs:ize)) &
+ grad2par_x_n(iy,ix) *f_x(iy,ix,izs:ize) &
+ grad2par_y_n(iy,ix) *f_y(iy,ix,izs:ize)
end do
end do
end subroutine grad2par_n2n_fd4
!-----------------------------------------------------------------------------
!-------------------------------------------------------------------------------
subroutine grad2par_v2v_fd4(f,flu,frd,f_lapl)
! flu and frd are useless here bum passing them anyway to respect the code structure in numerics
use prec_const
IMPLICIT none
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), INTENT(in) :: f
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2), intent(in) :: flu, frd
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_lapl
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_zz, f_yy, f_xx
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_zy, f_yz, f_zx, f_xz, f_xy, f_yx
real(dp), DIMENSION(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_x, f_y, f_z
integer :: ix,iy
call gradxx_n2n_fd4(f ,f_xx)
call gradyy_n2n_fd4(f ,f_yy)
call gradzz_n2n_fd4(f ,f_zz)
call gradx_n2n_fd4 (f ,f_x )
call grady_n2n_fd4 (f ,f_y )
call gradz_n2n_fd4 (f ,f_z )
call grady_n2n_fd4 (f_z,f_yz)
call gradz_n2n_fd4 (f_y,f_zy)
call gradx_n2n_fd4 (f_z,f_xz)
call gradz_n2n_fd4 (f_x,f_zx)
call grady_n2n_fd4 (f_x,f_yx)
call gradx_n2n_fd4 (f_y,f_xy)
f_lapl(:,:,:) = nan_
do ix = ixs,ixe
do iy = iys, iye
f_lapl(iy,ix,izs:ize) = f_zz(iy,ix,izs:ize) &
+ grad2par_yy_v(iy,ix) *f_yy(iy,ix,izs:ize) &
+ grad2par_xx_v(iy,ix) *f_xx(iy,ix,izs:ize) &
+ grad2par_zy_v(iy,ix) *(f_zy(iy,ix,izs:ize)+f_yz(iy,ix,izs:ize)) &
+ grad2par_zx_v(iy,ix) *(f_zx(iy,ix,izs:ize)+f_xz(iy,ix,izs:ize)) &
+ grad2par_xy_v(iy,ix) *(f_xy(iy,ix,izs:ize)+f_yx(iy,ix,izs:ize)) &
+ grad2par_x_v(iy,ix) *f_x(iy,ix,izs:ize) &
+ grad2par_y_v(iy,ix) *f_y(iy,ix,izs:ize)
end do
end do
end subroutine grad2par_v2v_fd4
!-----------------------------------------------------------------------------
!*******************************************
!************** INTERPOLATION **************
!*******************************************
subroutine interp_v2n_fd4 ( f, flu, frd , f_int_v2n )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg),intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_int_v2n
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2), intent(in) :: flu, frd
integer :: iy, iz
f_int_v2n(:,:,:) = nan_
do iy=iys,iye
do iz=izs,ize
f_int_v2n(iy,:,iz) = coef_int(1)*(coef_int(1)*f(iy-1, :, iz-1) &
+ coef_int(2)*f(iy-1, : , iz) &
+ coef_int(3)*f(iy-1, :, iz +1) &
+ coef_int(4)*f(iy-1, :, iz+2) )&
+ coef_int(2)*( coef_int(1)*f(iy, :, iz-1) &
+ coef_int(2)*f(iy, :, iz) &
+ coef_int(3)*f(iy, :, iz+1) &
+ coef_int(4)*f(iy, :, iz+2) )&
+ coef_int(3)*( coef_int(1)*f(iy+1, :, iz-1) &
+ coef_int(2)*f(iy+1, :, iz) &
+ coef_int(3)*f(iy+1, :, iz+1) &
+ coef_int(4)*f(iy+1, :, iz+2) )&
+ coef_int(4)*( coef_int(1)*f(iy+2, :, iz-1) &
+ coef_int(2)*f(iy+2, :, iz) &
+ coef_int(3)*f(iy+2, :, iz+1) &
+ coef_int(4)*f(iy+2, :, iz+2) )
end do
end do
end subroutine interp_v2n_fd4
!--------------------------------------------------------------------------------
subroutine interp_n2v_fd4( f, flu, frd , f_int_n2v )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg),intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_int_n2v
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2), intent(in) :: flu, frd
integer :: iy, iz
f_int_n2v(:,:,:) = nan_
do iy=iys,iye
do iz=izs,ize
f_int_n2v(iy,:,iz) = coef_int(1)*( coef_int(1)*f(iy-2, :, iz-2) &
+ coef_int(2)*f(iy-2, :, iz-1) &
+ coef_int(3)*f(iy-2, :, iz ) &
+ coef_int(4)*f(iy-2, :, iz+1) )&
+ coef_int(2)*( coef_int(1)*f(iy-1, :, iz-2) &
+ coef_int(2)*f(iy-1, :, iz-1) &
+ coef_int(3)*f(iy-1, :, iz) &
+ coef_int(4)*f(iy-1, :, iz+1) )&
+ coef_int(3)*( coef_int(1)*f(iy, :, iz-2) &
+ coef_int(2)*f(iy, :, iz-1) &
+ coef_int(3)*f(iy, :, iz) &
+ coef_int(4)*f(iy, :, iz+1) )&
+coef_int(4)*( coef_int(1)*f(iy+1, :, iz-2) &
+ coef_int(2)*f(iy+1, :, iz-1) &
+ coef_int(3)*f(iy+1, :, iz) &
+ coef_int(4)*f(iy+1, :, iz+1) )
end do
end do
end subroutine interp_n2v_fd4
!-----------------------------------------------------------------------------
!********************************************************
!**************** CURVATURE OPERATOR ********************
!********************************************************
!NB : the curvature operator here should be C/B not just C
!in reality since Im neglecting B everywhere it doesnt really matters
!!!! 3 curvature operators - n2n, n2v, v2n (in y) - NB cur_x_v is actually useless
!Curvature operator in s-alpha field representation
subroutine curv_n2n_fd4( f , C_f )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C_f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_x, f_y
integer :: iy, ix
call gradx_n2n_fd4( f, f_x)
call grady_n2n_fd4( f, f_y)
C_f(:,:,:) = nan_
do ix=ixs,ixe
do iy=iys,iye
C_f(iy,ix,izs:ize) = f_x(iy,ix,izs:ize)*cur_x_n(iy) + f_y(iy,ix,izs:ize)*cur_y_n(iy,ix)
end do
end do
end subroutine curv_n2n_fd4
subroutine curv_n2n_rec(f, C_f)
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C_f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_y
integer :: iy, ix
call grady_n2n_fd4( f, f_y)
C_f(:,:,:) = nan_
do ix=ixs,ixe
do iy=iys,iye
C_f(iy,ix,izs:ize) = f_y(iy,ix,izs:ize)*cur_y_n(iy,ix)
end do
end do
end subroutine curv_n2n_rec
!-------------------------------------------------------------
subroutine curv_curv_n2n_fd4 (f , f_cc)
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: f_cc
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_xx, f_yy, f_x, f_xy
integer :: iy, ix
call gradxx_n2n_fd4(f,f_xx)
call gradyy_n2n_fd4(f,f_yy)
call gradx_n2n_fd4 (f ,f_x )
call grady_n2n_fd4 (f_x,f_xy)
f_cc(:,:,:) = nan_
do ix=ixs,ixe
do iy=iys,iye
f_cc(iy,ix,izs:ize) = f_xx(iy,ix,izs:ize)*ccur_xx_n(iy) &
+ f_yy(iy,ix,izs:ize)*ccur_yy_n(iy,ix) &
+ f_xy(iy,ix,izs:ize)*ccur_xy_n(iy,ix)
end do
end do
end subroutine curv_curv_n2n_fd4
subroutine curv_curv_n2n_rec (f , f_cc)
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: f_cc
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_yy
integer :: iy, ix
call gradyy_n2n_fd4(f,f_yy)
f_cc(:,:,:) = nan_
do ix=ixs,ixe
do iy=iys,iye
f_cc(iy,ix,izs:ize) = f_yy(iy,ix,izs:ize)*ccur_yy_n(iy,ix)
end do
end do
end subroutine curv_curv_n2n_rec
!--------------------------------------------------------------
subroutine curv_v2v_fd4( f , C_f )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: C_f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_x, f_y
integer :: iy, ix
call gradx_n2n_fd4( f, f_x)
call grady_n2n_fd4( f, f_y)
C_f(:,:,:) = nan_
do ix=ixs,ixe
do iy=iys,iye
C_f(iy,ix,izs:ize) = f_x(iy,ix,izs:ize)*cur_x_v(iy) + f_y(iy,ix,izs:ize)*cur_y_v(iy,ix)
end do
end do
end subroutine curv_v2v_fd4
subroutine curv_v2v_rec( f , C_f )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: C_f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_y
integer :: iy, ix
call grady_n2n_fd4( f, f_y)
C_f(:,:,:) = nan_
do ix=ixs,ixe
do iy=iys,iye
C_f(iy,ix,izs:ize) = f_y(iy,ix,izs:ize)*cur_y_v(iy,ix)
end do
end do
end subroutine curv_v2v_rec
!----------------------------------------------------------------
subroutine curv_n2v_fd4( f , C_f )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C_f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_x, f_y, f_x_v
integer :: iy, ix
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
call gradx_n2n_fd4( f, f_x)
call grady_n2v_fd4( f, f_y)
call interp_n2v_fd4( f_x, buff, buff , f_x_v)
C_f(:,:,:) = nan_
do ix=ixs,ixe
do iy=iys,iye
C_f(iy,ix,izs:ize) = f_x_v(iy,ix,izs:ize)*cur_x_v(iy) + f_y(iy,ix,izs:ize)*cur_y_v(iy,ix)
end do
end do
end subroutine curv_n2v_fd4
subroutine curv_n2v_rec( f , C_f )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C_f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_y
integer :: iy, ix
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
call grady_n2v_fd4( f, f_y)
C_f(:,:,:) = nan_
do ix=ixs,ixe
do iy=iys,iye
C_f(iy,ix,izs:ize) = f_y(iy,ix,izs:ize)*cur_y_v(iy,ix)
end do
end do
end subroutine curv_n2v_rec
!-----------------------------------------------------------------
subroutine curv_v2n_fd4( f , C_f )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C_f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_x, f_y, f_x_n
integer :: iy, ix
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
call gradx_n2n_fd4( f, f_x)
call interp_v2n_fd4(f_x,buff,buff,f_x_n)
call grady_v2n_fd4( f,f_y )
C_f (:,:,:) =nan_
do ix=ixs,ixe
do iy=iys,iye
C_f(iy,ix,izs:ize) = f_x_n(iy,ix,izs:ize)*cur_x_n(iy) + f_y(iy,ix,izs:ize)*cur_y_n(iy,ix)
end do
end do
end subroutine curv_v2n_fd4
subroutine curv_v2n_rec( f , C_f )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C_f
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_y
integer :: iy, ix
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
call grady_v2n_fd4( f,f_y )
C_f (:,:,:) =nan_
do ix=ixs,ixe
do iy=iys,iye
C_f(iy,ix,izs:ize) = f_y(iy,ix,izs:ize)*cur_y_n(iy,ix)
end do
end do
end subroutine curv_v2n_rec
!-----------------------------------------------------------------
!*************************************************
!************* POISSON BRACKET *******************
!*************************************************
subroutine arakw_h(A,B, C)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: A, B
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: C
real(dp):: z
INTEGER m,n,k
z = - 1._dp/(12._dp*deltax*deltay)
C(:,:,:) = nan_
DO k=izs,ize
DO n=ixs,ixe
DO m=iys,iye
C(m,n,k) = z*( &
(A(m ,n-1,k) + A(m+1,n-1,k) - A(m ,n+1,k) - A(m+1,n+1,k))*(B(m+1,n ,k) - B(m ,n ,k)) &
+(A(m-1,n-1,k) + A(m ,n-1,k) - A(m-1,n+1,k) - A(m ,n+1,k))*(B(m ,n ,k) - B(m-1,n ,k)) &
+(A(m+1,n ,k) + A(m+1,n+1,k) - A(m-1,n ,k) - A(m-1,n+1,k))*(B(m,n+1 ,k) - B(m ,n ,k)) &
+(A(m+1,n-1,k) + A(m+1,n ,k) - A(m-1,n-1,k) - A(m-1,n ,k))*(B(m,n ,k) - B(m ,n-1,k)) &
+(A(m+1,n ,k) - A(m ,n+1,k)) *(B(m+1,n+1,k) - B(m ,n ,k)) &
+(A(m ,n-1,k) - A(m-1,n ,k)) *(B(m ,n ,k) - B(m-1,n-1,k)) &
+(A(m ,n+1,k) - A(m-1,n ,k)) *(B(m-1,n+1,k) - B(m ,n ,k)) &
+(A(m+1,n ,k) - A(m ,n-1,k)) *(B(m ,n ,k) - B(m+1,n-1,k)) )
END DO
END DO
END DO
end subroutine arakw_h
!-----------------------------------------------------------------
subroutine arakw_2h(A,B, C)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: A, B
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: C
real(dp):: z
INTEGER m,n,k
z = - 1._dp/(12._dp*(2._dp*deltax)*(2._dp*deltay))
C(:,:,:) = nan_
DO k=izs,ize
DO n=ixs,ixe
DO m=iys,iye
C(m,n,k) = z*( &
(A(m ,n-2,k) + A(m+2,n-2,k) - A(m ,n+2,k) - A(m+2,n+2,k))*(B(m+2,n ,k) - B(m ,n ,k)) &
+(A(m-2,n-2,k) + A(m ,n-2,k) - A(m-2,n+2,k) - A(m ,n+2,k))*(B(m ,n ,k) - B(m-2,n ,k)) &
+(A(m+2,n ,k) + A(m+2,n+2,k) - A(m-2,n ,k) - A(m-2,n+2,k))*(B(m,n+2 ,k) - B(m ,n ,k)) &
+(A(m+2,n-2,k) + A(m+2,n ,k) - A(m-2,n-2,k) - A(m-2,n ,k))*(B(m,n ,k) - B(m ,n-2,k)) &
+(A(m+2,n ,k) - A(m ,n+2,k)) *(B(m+2,n+2,k) - B(m ,n ,k)) &
+(A(m ,n-2,k) - A(m-2,n ,k)) *(B(m ,n ,k) - B(m-2,n-2,k)) &
+(A(m ,n+2,k) - A(m-2,n ,k)) *(B(m-2,n+2,k) - B(m ,n ,k)) &
+(A(m+2,n ,k) - A(m ,n-2,k)) *(B(m ,n ,k) - B(m+2,n-2,k)) )
END DO
END DO
END DO
end subroutine arakw_2h
!-----------------------------------------------------------------
subroutine arakw_n2n_fd4(A,B, C)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: A, B
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: C
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg):: C_h, C_2h
INTEGER ix, iy, iz
call arakw_h(A,B, C_h)
call arakw_2h(A,B, C_2h)
C(:,:,:) = nan_
do iz=izs,ize
do ix=ixs,ixe
do iy = iys, iye
C(iy,ix,iz) = (4._dp/3._dp*C_h(iy,ix,iz) - 1._dp/3._dp*C_2h(iy,ix,iz))
end do
end do
end do
end subroutine arakw_n2n_fd4
!-----------------------------------------------------------------
subroutine arakw_n2v_fd4(A,B,C)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: A, B
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg):: A_x_n, A_x_v, A_y_v, B_x_v, B_y_v
INTEGER ix, iy, iz
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
! A is tipically the stream subroutine which defined on n grid and B are velocity defined on v grid
call gradx_n2n_fd4(A,A_x_n)
call interp_n2v_fd4(A_x_n,buff, buff,A_x_v)
call grady_n2v_fd4(A, A_y_v)
call gradx_n2n_fd4(B,B_x_v)
call grady_n2n_fd4(B,B_y_v)
C(:,:,:) = nan_
do iz=izs,ize
do ix=ixs,ixe
do iy = iys, iye
C(iy,ix,iz) = ( A_x_v(iy,ix,iz)*B_y_v(iy,ix,iz) - A_y_v(iy,ix,iz)*B_x_v(iy,ix,iz))
end do
end do
end do
end subroutine arakw_n2v_fd4
subroutine arakw_n2v_rec(A,B,C)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: A, B
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg):: B_x_n, B_x_v, B_y_v, A_x_v, A_y_v
INTEGER ix, iy, iz
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
! B is tipically the stream subroutine which is defined on n grid and A are velocities defined on v grid
call gradx_n2n_fd4(B,B_x_n)
call interp_n2v_fd4(B_x_n,buff, buff,B_x_v)
call grady_n2v_fd4(B, B_y_v)
call gradx_n2n_fd4(A,A_x_v)
call grady_n2n_fd4(A,A_y_v)
C(:,:,:) = nan_
do iz=izs,ize
do ix=ixs,ixe
do iy = iys, iye
C(iy,ix,iz) = ( A_x_v(iy,ix,iz)*B_y_v(iy,ix,iz) - A_y_v(iy,ix,iz)*B_x_v(iy,ix,iz))
end do
end do
end do
end subroutine arakw_n2v_rec
subroutine arakw_n2v_v2v_rec(A,B,C)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: A, B
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg):: A_x_n, A_x_v, A_y_v, B_x_v, B_y_v
INTEGER ix, iy, iz
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
! A: n -> v, B: v -> v
call gradx_n2n_fd4(A,A_x_n)
call interp_n2v_fd4(A_x_n,buff, buff,A_x_v)
call grady_n2v_fd4(A, A_y_v)
call gradx_n2n_fd4(B,B_x_v)
call grady_n2n_fd4(B,B_y_v)
C(:,:,:) = nan_
do iz=izs,ize
do ix=ixs,ixe
do iy = iys, iye
C(iy,ix,iz) = ( A_x_v(iy,ix,iz)*B_y_v(iy,ix,iz) - A_y_v(iy,ix,iz)*B_x_v(iy,ix,iz))
end do
end do
end do
end subroutine arakw_n2v_v2v_rec
subroutine arakw_n2n_v2n_rec(A,B,C)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: A, B
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg):: A_x_n, A_y_n, B_x_v, B_x_n, B_y_n
INTEGER ix, iy, iz
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
! A: n -> n, B: v -> n
call gradx_n2n_fd4(B,B_x_v)
call interp_v2n_fd4(B_x_v,buff, buff,B_x_n)
call grady_v2n_fd4(B, B_y_n)
call gradx_n2n_fd4(A,A_x_n)
call grady_n2n_fd4(A,A_y_n)
C(:,:,:) = nan_
do iz=izs,ize
do ix=ixs,ixe
do iy = iys, iye
C(iy,ix,iz) = ( A_x_n(iy,ix,iz)*B_y_n(iy,ix,iz) - A_y_n(iy,ix,iz)*B_x_n(iy,ix,iz))
end do
end do
end do
end subroutine arakw_n2n_v2n_rec
subroutine arakw_v2n_v2n_rec(A,B,C)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: A, B
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: C
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg):: A_x_v, A_x_n, A_y_n, B_x_v, B_x_n, B_y_n
INTEGER ix, iy, iz
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
! A: v -> n, B: v -> n
call gradx_n2n_fd4(A,A_x_v)
call interp_v2n_fd4(A_x_v,buff, buff,A_x_n)
call grady_v2n_fd4(A, A_y_n)
call gradx_n2n_fd4(B,B_x_v)
call interp_v2n_fd4(B_x_v,buff, buff,B_x_n)
call grady_v2n_fd4(B, B_y_n)
C(:,:,:) = nan_
do iz=izs,ize
do ix=ixs,ixe
do iy = iys, iye
C(iy,ix,iz) = ( A_x_n(iy,ix,iz)*B_y_n(iy,ix,iz) - A_y_n(iy,ix,iz)*B_x_n(iy,ix,iz))
end do
end do
end do
end subroutine arakw_v2n_v2n_rec
!-----------------------------------------------------------------
!******************************************************
!*************** GRAD CURV giroviscous ****************
!******************************************************
!only the one from v to v is needed because we have usually grad_par(C(strmf)+C(T)) in the v eq
! the strmf has to be first interpolated on the v grid and its BC updated
!THIS ROUTINES NEED TO BE CHANGED!!!!
subroutine gradpar_curv_n2v_fd4 (f_n, f_parC)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: f_n
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_parC
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_xx_n, f_xx, f_yy, f_xy, f_zx, f_zy, f_x_n, f_x, f_y
INTEGER ix, iy
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
! I dont know if this works
call gradxx_n2n_fd4(f_n,f_xx_n )
call interp_n2v_fd4(f_xx_n,buff,buff, f_xx)
call gradyy_n2v_fd4(f_n, f_yy)
call gradx_n2n_fd4(f_n,f_x_n)
call interp_n2v_fd4(f_x_n, buff, buff,f_x)
call grady_n2v_fd4(f_n,f_y)
call gradx_n2n_fd4(f_y,f_xy)
! it seems that numerically with FD doing dzy or dyz is the same
call gradyz_n2v_fd4(f_n, f_zy )
call gradz_n2v_fd4(f_x_n, f_zx)
f_parC(:,:,:) = nan_
do ix=ixs,ixe
do iy = iys, iye
f_parC(iy,ix,izs:ize) = gradcur_xx_v(iy,ix)*f_xx(iy,ix,izs:ize) &
+ gradcur_yy_v(iy,ix)*f_yy(iy,ix,izs:ize) &
+ (gradcur_xy_v(iy,ix) + gradcur_yx_v(iy,ix))*f_xy(iy,ix,izs:ize) &
+ gradcur_zx_v(iy,ix)*f_zx(iy,ix,izs:ize) &
+ gradcur_zy_v(iy,ix)*f_zy(iy,ix,izs:ize) &
+ gradcur_x_v (iy,ix)*f_x (iy,ix,izs:ize) &
+ gradcur_y_v (iy,ix)*f_y (iy,ix,izs:ize)
end do
end do
end subroutine gradpar_curv_n2v_fd4
subroutine gradpar_curv_n2v_rec (f_n, f_parC)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: f_n
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out) :: f_parC
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg) :: f_yy, f_xy, f_zy, f_y
INTEGER ix, iy
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
call gradyy_n2v_fd4(f_n, f_yy)
call grady_n2v_fd4(f_n,f_y)
call gradx_n2n_fd4(f_y,f_xy)
call gradyz_n2v_fd4(f_n, f_zy )
f_parC(:,:,:) = nan_
do ix=ixs,ixe
do iy = iys, iye
f_parC(iy,ix,izs:ize) = gradcur_yy_v(iy,ix)*f_yy(iy,ix,izs:ize) + gradcur_xy_v(iy,ix)*f_xy(iy,ix,izs:ize) + gradcur_zy_v(iy,ix)*f_zy(iy,ix,izs:ize)
end do
end do
end subroutine gradpar_curv_n2v_rec
!-----------------------------------------------------------------
! this operator is used in the expression for omega and it is applied to v. Therefore it goes from v to n
subroutine curv_gradpar_v2n_fd4(f_v, f_Cpar)
use prec_const
IMPLICIT NONE
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in):: f_v
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(out):: f_Cpar
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg):: f_xx_v, f_xx, f_yy, f_xy, f_zx, f_zy, f_x_v, f_x, f_y
INTEGER ix, iy
real(dp), dimension(iylg:ny_zg,ixsg:ixeg,2) :: buff
buff(:,:,:) = 0._dp
! I dont know if this works
call gradxx_n2n_fd4(f_v, f_xx_v)
call interp_v2n_fd4(f_xx_v,buff,buff, f_xx)
call gradyy_v2n_fd4(f_v,f_yy )
call gradx_n2n_fd4 (f_v,f_x_v)
call interp_v2n_fd4(f_x_v, buff, buff,f_x)
call grady_v2n_fd4 (f_v, f_y)
call gradx_n2n_fd4(f_y, f_xy)
! it seems that numerically with FD doing dzy or dyz is the same
call gradyz_v2n_fd4(f_v,f_zy)
call gradz_v2n_fd4(f_x_v,f_zx)
f_Cpar(:,:,:) = nan_
do ix=ixs,ixe
do iy = iys, iye
f_Cpar (iy,ix,izs:ize) = curgrad_xx_n(iy,ix)*f_xx(iy,ix,izs:ize) &
+ curgrad_yy_n(iy,ix)*f_yy(iy,ix,izs:ize) &
+ (curgrad_xy_n(iy,ix) + curgrad_yx_n(iy,ix))*f_xy(iy,ix,izs:ize) &
+ curgrad_zx_n(iy,ix)*f_zx(iy,ix,izs:ize) &
+ curgrad_zy_n(iy,ix)*f_zy(iy,ix,izs:ize) &
+ curgrad_x_n (iy,ix)*f_x (iy,ix,izs:ize) &
+ curgrad_y_n (iy,ix)*f_y (iy,ix,izs:ize)
end do
end do
end subroutine curv_gradpar_v2n_fd4
!******************************************************
!*************** Diffusion operators ****************
!******************************************************
SUBROUTINE scalxx_fa( f, diff , fxx )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f, diff
real(dp), dimension(iys:iye ,ixs:ixe ,izs:ize) , intent(out):: fxx
integer :: iy, ix, iz
!$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
do iz=izs,ize
do ix=ixs,ixe
do iy=iys,iye
fxx(iy,ix,iz) = deltaxisq * ( &
f(iy,ix-1,iz) * 0.5_DP * ( diff(iy,ix,iz) + diff(iy,ix-1,iz) ) &
-2.0_DP*f(iy,ix ,iz) * 0.5_DP * ( diff(iy,ix,iz) + 0.5_DP*( diff(iy,ix-1,iz) + diff(iy,ix+1,iz) ) ) &
+ f(iy,ix+1,iz) * 0.5_DP * ( diff(iy,ix,iz) + diff(iy,ix+1,iz) ) &
)
end do
end do
end do
!$OMP END DO
END SUBROUTINE scalxx_fa
SUBROUTINE scalyy_fa( f, diff , fyy )
use prec_const
implicit none
real(dp), dimension(iysg:iyeg,ixsg:ixeg,izsg:izeg), intent(in) :: f, diff
real(dp), dimension(iys:iye ,ixs:ixe ,izs:ize) , intent(out):: fyy
integer :: iy, ix, iz
!$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
do iz=izs,ize
do ix=ixs,ixe
do iy=iys,iye
fyy(iy,ix,iz) = deltayisq * ( &
f(iy-1,ix,iz) * 0.5_DP * ( diff(iy,ix,iz) + diff(iy-1,ix,iz) ) &
-2.0_DP*f(iy ,ix,iz) * 0.5_DP * ( diff(iy,ix,iz) + 0.5_DP*( diff(iy-1,ix,iz) + diff(iy+1,ix,iz) ) ) &
+ f(iy+1,ix,iz) * 0.5_DP * ( diff(iy,ix,iz) + diff(iy+1,ix,iz) ) &
)
end do
end do
end do
!$OMP END DO
END SUBROUTINE scalyy_fa
c4science · Help