!= Module xyz_deriv_c4_module ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: xyz_deriv_c4_module.f90,v 1.4 2014-03-04 05:55:07 sugiyama Exp $ ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] module xyz_deriv_c4_module != 3 次元 (xyz 方向) 不等間隔交互格子 微分演算モジュール 4 次精度版 ! ! xyz_deriv_module は, 3 次元 (xyz 方向) 不等間隔交互格子を用いた有限差分 ! 法に基づく数値モデルのための, 微分演算を行う Fortran90 関数を提供する. ! 微分演算は 4 次精度中心差分を用いて行う. ! ! 差分式は格子間隔を一定とした ! ! du/dx(i) = 9/8*[u(i)-u(i-1)]/dx(i) - 1/24*[u(i+1)-u(i-2)]/dx(i) ! ! を用いている. use dc_types, only : DP use dc_message, only : MessageNotify use gridset, only : imin, imax, jmin, jmax, kmin, kmax use axesset, only : x_dx, p_dx, y_dy, q_dy, z_dz, r_dz implicit none private public :: xyz_c4dx_pyz, pyz_c4dx_xyz, pyr_c4dx_xyr, pqz_c4dx_xqz, & & xyz_c4dy_xqz, xqz_c4dy_xyz, pqz_c4dy_pyz, xqr_c4dy_xyr, & & xyz_c4dz_xyr, xyr_c4dz_xyz, pyr_c4dz_pyz, xqr_c4dz_xqz interface pyz_c4dx_xyz module procedure paa_c4dx_xaa end interface interface pyr_c4dx_xyr module procedure paa_c4dx_xaa end interface interface pqz_c4dx_xqz module procedure paa_c4dx_xaa end interface interface xqz_c4dy_xyz module procedure aqa_c4dy_aya end interface interface pqz_c4dy_pyz module procedure aqa_c4dy_aya end interface interface xqr_c4dy_xyr module procedure aqa_c4dy_aya end interface interface xyr_c4dz_xyz module procedure aar_c4dz_aaz end interface interface pyr_c4dz_pyz module procedure aar_c4dz_aaz end interface interface xqr_c4dz_xqz module procedure aar_c4dz_aaz end interface contains !-------------------------------------------------------------------- function xyz_c4dx_pyz(pyz_Var) ! 半整数格子点上の 1 階微分を計算する real(DP),intent(in) :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_c4dx_pyz(imin:imax,jmin:jmax,kmin:kmax) integer :: ix ! 1 階微分の計算 ! do ix = imin+2, imax-1 xyz_c4dx_pyz(ix,:,:) = & & ((pyz_Var(ix,:,:) - pyz_Var(ix-1,:,:))*9.0d0/(8.0d0*x_dx(ix)) & & - ((pyz_Var(ix+1,:,:) - pyz_Var(ix-2,:,:))/(24.0d0*x_dx(ix)))) end do ! 値が決まらない格子点の値を適当に決める ! 計算領域に滲み出したらすぐにわかるように, 大きな値を入れておく ! xyz_c4dx_pyz(imin:imin+1,:,:) = 1.0d10 xyz_c4dx_pyz(imax,:,:) = 1.0d10 end function xyz_c4dx_pyz !-------------------------------------------------------------------- function paa_c4dx_xaa(xaa_Var) ! 半整数格子点上の 1 階微分を計算する real(DP),intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: paa_c4dx_xaa(imin:imax,jmin:jmax,kmin:kmax) integer :: ix ! 1 階微分の計算 ! do ix = imin+1, imax-2 paa_c4dx_xaa(ix,:,:) = & & ((xaa_Var(ix+1,:,:) - xaa_Var(ix,:,:))*9.0d0/(8.0d0*p_dx(ix)) & & - ((xaa_Var(ix+2,:,:) - xaa_Var(ix-1,:,:))/(24.0d0*p_dx(ix)))) end do ! 値が決まらない格子点の値を適当に決める ! 計算領域に滲み出したらすぐにわかるように, 大きな値を入れておく ! paa_c4dx_xaa(imin,:,:) = 1.0d10 paa_c4dx_xaa(imax-1:imax,:,:) = 1.0d10 end function paa_c4dx_xaa !-------------------------------------------------------------------- function xyz_c4dy_xqz(xqz_Var) ! 半整数格子点上の 1 階微分を計算する real(DP),intent(in) :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_c4dy_xqz(imin:imax,jmin:jmax,kmin:kmax) integer :: jy ! y 方向一様な場合は微分はゼロ. if (jmin == jmax) then xyz_c4dy_xqz = 0.0d0 return end if ! 1 階微分の計算 ! do jy = jmin+2, jmax-1 xyz_c4dy_xqz(:,jy,:) = & & ((xqz_Var(:,jy,:) - xqz_Var(:,jy-1,:))*9.0d0/(8.0d0*y_dy(jy)) & & - ((xqz_Var(:,jy+1,:) - xqz_Var(:,jy-2,:))/(24.0d0*y_dy(jy)))) end do ! 値が決まらない格子点の値を適当に決める ! 計算領域に滲み出したらすぐにわかるように, 大きな値を入れておく ! xyz_c4dy_xqz(:,jmin:jmin+1,:) = 1.0d10 xyz_c4dy_xqz(:,jmax,:) = 1.0d10 end function xyz_c4dy_xqz !-------------------------------------------------------------------- function aqa_c4dy_aya(aya_Var) ! 半整数格子点上の 1 階微分を計算する real(DP),intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: aqa_c4dy_aya(imin:imax,jmin:jmax,kmin:kmax) integer :: jy ! y 方向一様な場合は微分はゼロ. if (jmin == jmax) then aqa_c4dy_aya = 0.0d0 return end if ! 1 階微分の計算 ! do jy = jmin+1, jmax-2 aqa_c4dy_aya(:,jy,:) = & & ((aya_Var(:,jy+1,:) - aya_Var(:,jy,:))*9.0d0/(8.0d0*q_dy(jy)) & & - ((aya_Var(:,jy+2,:) - aya_Var(:,jy-1,:))/(24.0d0*q_dy(jy)))) end do ! 値が決まらない格子点の値を適当に決める ! 計算領域に滲み出したらすぐにわかるように, 大きな値を入れておく ! aqa_c4dy_aya(:,jmin,:) = 1.0d10 aqa_c4dy_aya(:,jmax-1:jmax,:) = 1.0d10 end function aqa_c4dy_aya !-------------------------------------------------------------------- function xyz_c4dz_xyr(xyr_Var) ! 半整数格子点上の 1 階微分を計算する real(DP),intent(in) :: xyr_Var(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_c4dz_xyr(imin:imax,jmin:jmax,kmin:kmax) integer :: kz ! 1 階微分の計算 ! do kz = kmin+2, kmax-1 xyz_c4dz_xyr(:,:,kz) = & & ((xyr_Var(:,:,kz) - xyr_Var(:,:,kz-1))*9.0d0/(8.0d0*z_dz(kz)) & & - ((xyr_Var(:,:,kz+1) - xyr_Var(:,:,kz-2))/(24.0d0*z_dz(kz)))) end do ! 値が決まらない格子点の値を適当に決める ! 計算領域に滲み出したらすぐにわかるように, 大きな値を入れておく ! xyz_c4dz_xyr(:,:,kmin:kmin+1) = 1.0d10 xyz_c4dz_xyr(:,:,kmax) = 1.0d10 end function xyz_c4dz_xyr !-------------------------------------------------------------------- function aar_c4dz_aaz(aaz_Var) ! 半整数格子点上の 1 階微分を計算する real(DP),intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: aar_c4dz_aaz(imin:imax,jmin:jmax,kmin:kmax) integer :: kz ! 1 階微分の計算 ! * r_c4dz_z を用いて計算 ! do kz = kmin+1, kmax-2 aar_c4dz_aaz(:,:,kz) = & & ((aaz_Var(:,:,kz+1) - aaz_Var(:,:,kz))*9.0d0/(8.0d0*r_dz(kz)) & & - ((aaz_Var(:,:,kz+2) - aaz_Var(:,:,kz-1))/(24.0d0*r_dz(kz)))) end do ! 値が決まらない格子点の値を適当に決める ! 計算領域に滲み出したらすぐにわかるように, 大きな値を入れておく ! aar_c4dz_aaz(:,:,kmin) = 1.0d10 aar_c4dz_aaz(:,:,kmax-1:kmax) = 1.0d10 end function aar_c4dz_aaz !-------------------------------------------------------------------- end module xyz_deriv_c4_module