[cig-commits] [commit] Hiro_latest: Use structure for insulated boundary (84b3344)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Nov 18 16:21:53 PST 2013
Repository : ssh://geoshell/calypso
On branch : Hiro_latest
Link : https://github.com/geodynamics/calypso/compare/93e9f8f974c7a247c8f02e54ec18de063f86c8fb...3c548304673360ddedd7d68c8095b3fb74a2b9ce
>---------------------------------------------------------------
commit 84b334495d7b3ff1b5259c383ec6b026725dbfca
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Tue Nov 12 22:00:45 2013 -0800
Use structure for insulated boundary
>---------------------------------------------------------------
84b334495d7b3ff1b5259c383ec6b026725dbfca
.../MHD_src/sph_MHD/cal_sol_sph_fluid_crank.f90 | 6 +-
.../MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90 | 72 ----------
.../MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90 | 154 ++++++++++++---------
.../MHD_src/sph_MHD/cal_sph_exp_nod_icb_ins.f90 | 146 +++++++++++--------
.../MHD_src/sph_MHD/const_sph_diffusion.f90 | 10 +-
.../MHD_src/sph_MHD/const_sph_radial_grad.f90 | 14 +-
.../MHD_src/sph_MHD/const_sph_rotation.f90 | 24 ++--
.../MHD_src/sph_MHD/m_addresses_trans_sph_MHD.f90 | 4 +-
.../MHD_src/sph_MHD/set_bc_sph_mhd.f90 | 68 +++++++--
.../MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90 | 4 +-
10 files changed, 267 insertions(+), 235 deletions(-)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sol_sph_fluid_crank.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sol_sph_fluid_crank.f90
index ee71a05..04f41ab 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sol_sph_fluid_crank.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sol_sph_fluid_crank.f90
@@ -195,7 +195,8 @@
call delete_zero_degree_comp(ipol%i_magne)
!
if(sph_bc_B%iflag_icb .eq. iflag_sph_insulator) then
- call cal_sph_nod_icb_ins_mag2(ipol%i_magne)
+ call cal_sph_nod_icb_ins_mag2(nidx_rj(2), sph_bc_B%kr_in, &
+ & ipol%i_magne)
else if(sph_bc_B%iflag_icb .eq. iflag_radial_magne) then
call cal_sph_nod_icb_qvc_mag2(ipol%i_magne)
end if
@@ -203,7 +204,8 @@
if(sph_bc_B%iflag_cmb .eq. iflag_radial_magne) then
call cal_sph_nod_cmb_qvc_mag2(ipol%i_magne)
else
- call cal_sph_nod_cmb_ins_mag2(ipol%i_magne)
+ call cal_sph_nod_cmb_ins_mag2(nidx_rj(2), sph_bc_B%kr_out, &
+ & ipol%i_magne)
end if
!
call lubksb_3band_mul(np_smp, idx_rj_smp_stack(0,2), &
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90
deleted file mode 100644
index cf55107..0000000
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90
+++ /dev/null
@@ -1,72 +0,0 @@
-!> @file cal_sph_bc_fdm_matrix.f90
-!! module cal_sph_bc_fdm_matrix
-!!
-!! @author H. Matsui
-!! @date Written on May, 2003
-!
-!!> @brief calculate FDM matrices for boundaries
-!!
-!!@verbatim
-!! subroutine s_cal_sph_bc_fdm_matrices
-!!@endverbatim
-!
- module cal_sph_bc_fdm_matrix
-!
- use m_precision
-!
- use m_machine_parameter
-!
- implicit none
-!
-! -----------------------------------------------------------------------
-!
- contains
-!
-! -----------------------------------------------------------------------
-!
- subroutine s_cal_sph_bc_fdm_matrices
-!
- use m_spheric_parameter
- use m_coef_fdm_fixed_ICB
- use m_coef_fdm_fixed_CMB
- use m_coef_fdm_free_ICB
- use m_coef_fdm_free_CMB
- use m_coef_fdm_to_center
- use cal_fdm_coefs_4_boundaries
-!
-!
- call cal_fdm2_coef_fix_fld_ICB(radius_1d_rj_r(nlayer_ICB), &
- & coef_fdm_fix_ICB_2)
- call cal_fdm2_coef_fix_df_ICB(radius_1d_rj_r(nlayer_ICB), &
- & coef_fdm_fix_dr_ICB_2)
-!
- call cal_fdm2_coef_fix_fld_CMB(radius_1d_rj_r(nlayer_CMB-2), &
- & coef_fdm_fix_CMB_2)
- call cal_fdm2_coef_fix_df_CMB(radius_1d_rj_r(nlayer_CMB-1), &
- & coef_fdm_fix_dr_CMB_2)
-!
-!
- call cal_2nd_ICB_free_vp_bc_fdm(radius_1d_rj_r(nlayer_ICB))
- call cal_2nd_ICB_free_vt_bc_fdm(radius_1d_rj_r(nlayer_ICB))
-!
- call cal_2nd_CMB_free_vp_bc_fdm(radius_1d_rj_r(nlayer_CMB-1))
- call cal_2nd_CMB_free_vt_bc_fdm(radius_1d_rj_r(nlayer_CMB-1))
-!
- call cal_2nd_to_center_fixed_fdm(radius_1d_rj_r(1))
- call cal_2nd_to_center_fix_df_fdm(radius_1d_rj_r(1))
-!
- if (iflag_debug .eq. iflag_full_msg) then
- call check_coef_fdm_fix_dr_ICB
- call check_coef_fdm_fix_dr_CMB
- call check_coef_fdm_free_ICB
- call check_coef_fdm_free_CMB
- call check_coef_fdm_fix_dr_2ctr
- end if
-!
-! call cal_sph_bc_2nd_ele_fdm_mat
-!
- end subroutine s_cal_sph_bc_fdm_matrices
-!
-! -----------------------------------------------------------------------
-!
- end module cal_sph_bc_fdm_matrix
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90
index ae5a46f..31f1578 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90
@@ -3,20 +3,32 @@
!!
!!@author H. Matsui
!!@date Programmed in Jan., 2010
+!!@date Modified in Nov., 2013
!
!>@brief Set insulated magnetic boundary condition for CMB
!!
!!@verbatim
-!! subroutine cal_sph_nod_cmb_ins_b_and_j(is_fld, is_rot)
-!! subroutine cal_sph_nod_cmb_ins_mag2(is_fld)
+!! subroutine cal_sph_nod_cmb_ins_b_and_j(jmax, kr_out, &
+!! & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, is_fld, is_rot)
+!! subroutine cal_sph_nod_cmb_ins_mag2(jmax, kr_out, is_fld)
!!
-!! subroutine cal_sph_nod_cmb_ins_vp_rot2(is_fld, is_rot)
-!! subroutine cal_sph_nod_cmb_ins_rot2(is_fld, is_rot)
-!! subroutine cal_sph_nod_cmb_ins_diffuse2(coef_d, &
+!! subroutine cal_sph_nod_cmb_ins_vp_rot2(jmax, kr_out, &
+!! & is_fld, is_rot)
+!! subroutine cal_sph_nod_cmb_ins_rot2(jmax, kr_out, &
+!! & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, is_fld, is_rot)
+!! subroutine cal_sph_nod_cmb_ins_diffuse2(jmax, kr_out, &
+!! & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, coef_d, &
!! & is_fld, is_diffuse)
!!@endverbatim
!!
-!!@n @param coef_d Coefficient for diffusion term
+!!@n @param jmax Number of modes for spherical harmonics @f$L*(L+2)@f$
+!!@n @param kr_out Radial ID for outer boundary
+!!@n @param fdm2_fix_fld_CMB(0:2,3)
+!!! Matrix to evaluate radial derivative at CMB with fiexed field
+!!@n @param fdm2_fix_dr_CMB(-1:1,3)
+!! Matrix to evaluate field at CMB with fiexed radial derivative
+!!
+!!@n @param coef_d Coefficient for diffusion term
!!@n @param is_fld Field address of input field
!!@n @param is_rot Field address for curl of field
!!@n @param is_diffuse Field address for diffusion of field
@@ -39,38 +51,39 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_ins_b_and_j(is_fld, is_rot)
-!
- use m_coef_fdm_fixed_CMB
+ subroutine cal_sph_nod_cmb_ins_b_and_j(jmax, kr_out, &
+ & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, is_fld, is_rot)
!
- integer(kind = kint), intent(in) :: is_fld
- integer(kind = kint), intent(in) :: is_rot
+ integer(kind = kint), intent(in) :: jmax, kr_out
+ integer(kind = kint), intent(in) :: is_fld, is_rot
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_CMB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
!
real(kind = kreal) :: d1s_dr1,d2s_dr2, d1t_dr1
integer(kind = kint) :: j, inod, i_n1, i_n2
!
!
!$omp parallel do private(inod,i_n1,i_n2,j,d1s_dr1,d2s_dr2,d1t_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_CMB-1) * nidx_rj(2)
- i_n1 = inod - nidx_rj(2)
- i_n2 = i_n1 - nidx_rj(2)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ i_n1 = inod - jmax
+ i_n2 = i_n1 - jmax
!
- d1s_dr1 = - g_sph_rj(j,1) * ar_1d_rj(nlayer_CMB,1) &
+ d1s_dr1 = - g_sph_rj(j,1) * ar_1d_rj(kr_out,1) &
& * d_rj(inod,is_fld)
- d2s_dr2 = coef_fdm_fix_dr_CMB_2(-1,3) * d_rj(i_n1,is_fld ) &
- & + (coef_fdm_fix_dr_CMB_2( 0,3) &
- & - g_sph_rj(j,1)*ar_1d_rj(nlayer_CMB,1) &
- & *coef_fdm_fix_dr_CMB_2( 1,3)) * d_rj(inod,is_fld )
- d1t_dr1 = coef_fdm_fix_CMB_2(2,2) * d_rj(i_n2,is_fld+2) &
- & + coef_fdm_fix_CMB_2(1,2) * d_rj(i_n1,is_fld+2)
+ d2s_dr2 = fdm2_fix_dr_CMB(-1,3) * d_rj(i_n1,is_fld ) &
+ & + (fdm2_fix_dr_CMB( 0,3) &
+ & - g_sph_rj(j,1)*ar_1d_rj(kr_out,1) &
+ & *fdm2_fix_dr_CMB( 1,3)) * d_rj(inod,is_fld )
+ d1t_dr1 = fdm2_fix_fld_CMB(2,2) * d_rj(i_n2,is_fld+2) &
+ & + fdm2_fix_fld_CMB(1,2) * d_rj(i_n1,is_fld+2)
!
d_rj(inod,is_fld+1) = d1s_dr1
d_rj(inod,is_fld+2) = zero
d_rj(inod,is_rot ) = zero
d_rj(inod,is_rot+1) = d1t_dr1
d_rj(inod,is_rot+2) = - ( d2s_dr2 - g_sph_rj(j,3) &
- & * ar_1d_rj(nlayer_CMB,2)*d_rj(inod,is_fld ) )
+ & * ar_1d_rj(kr_out,2)*d_rj(inod,is_fld ) )
end do
!$omp end parallel do
!
@@ -78,8 +91,9 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_ins_mag2(is_fld)
+ subroutine cal_sph_nod_cmb_ins_mag2(jmax, kr_out, is_fld)
!
+ integer(kind = kint), intent(in) :: jmax, kr_out
integer(kind = kint), intent(in) :: is_fld
!
real(kind = kreal) :: d1s_dr1
@@ -87,9 +101,9 @@
!
!
!$omp parallel do private(inod,d1s_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_CMB-1) * nidx_rj(2)
- d1s_dr1 = - g_sph_rj(j,1) * ar_1d_rj(nlayer_CMB,1) &
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ d1s_dr1 = - g_sph_rj(j,1) * ar_1d_rj(kr_out,1) &
& * d_rj(inod,is_fld)
!
d_rj(inod,is_fld+1) = d1s_dr1
@@ -101,8 +115,10 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_ins_vp_rot2(is_fld, is_rot)
+ subroutine cal_sph_nod_cmb_ins_vp_rot2(jmax, kr_out, &
+ & is_fld, is_rot)
!
+ integer(kind = kint), intent(in) :: jmax, kr_out
integer(kind = kint), intent(in) :: is_fld, is_rot
!
real(kind = kreal) :: d1t_dr1
@@ -110,9 +126,9 @@
!
!
!$omp parallel do private(inod,d1t_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_CMB-1) * nidx_rj(2)
- d1t_dr1 = - g_sph_rj(j,1) * ar_1d_rj(nlayer_CMB,1) &
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ d1t_dr1 = - g_sph_rj(j,1) * ar_1d_rj(kr_out,1) &
& * d_rj(inod,is_fld+2)
!
d_rj(inod,is_rot ) = d_rj(inod,is_fld+2)
@@ -126,35 +142,36 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_ins_rot2(is_fld, is_rot)
-!
- use m_coef_fdm_fixed_CMB
+ subroutine cal_sph_nod_cmb_ins_rot2(jmax, kr_out, &
+ & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, is_fld, is_rot)
!
- integer(kind = kint), intent(in) :: is_fld
- integer(kind = kint), intent(in) :: is_rot
+ integer(kind = kint), intent(in) :: jmax, kr_out
+ integer(kind = kint), intent(in) :: is_fld, is_rot
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_CMB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
!
real(kind = kreal) :: d2s_dr2, d1t_dr1
integer(kind = kint) :: j, inod, i_n1, i_n2
!
!
!$omp parallel do private(inod,i_n1,i_n2,j,d2s_dr2,d1t_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_CMB-1) * nidx_rj(2)
- i_n1 = inod - nidx_rj(2)
- i_n2 = i_n1 - nidx_rj(2)
-!
- d2s_dr2 = coef_fdm_fix_dr_CMB_2(-1,3) * d_rj(i_n1,is_fld ) &
- & + (coef_fdm_fix_dr_CMB_2( 0,3) &
- & - g_sph_rj(j,1)*ar_1d_rj(nlayer_CMB,1) &
- & *coef_fdm_fix_dr_CMB_2( 1,3)) * d_rj(inod,is_fld )
- d1t_dr1 = coef_fdm_fix_CMB_2(2,2) * d_rj(i_n2,is_fld+2) &
- & + coef_fdm_fix_CMB_2(1,2) * d_rj(i_n1,is_fld+2) &
- & + coef_fdm_fix_CMB_2(0,2) * d_rj(inod,is_fld+2)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ i_n1 = inod - jmax
+ i_n2 = i_n1 - jmax
+!
+ d2s_dr2 = fdm2_fix_dr_CMB(-1,3) * d_rj(i_n1,is_fld ) &
+ & + (fdm2_fix_dr_CMB( 0,3) &
+ & - g_sph_rj(j,1)*ar_1d_rj(kr_out,1) &
+ & *fdm2_fix_dr_CMB( 1,3)) * d_rj(inod,is_fld )
+ d1t_dr1 = fdm2_fix_fld_CMB(2,2) * d_rj(i_n2,is_fld+2) &
+ & + fdm2_fix_fld_CMB(1,2) * d_rj(i_n1,is_fld+2) &
+ & + fdm2_fix_fld_CMB(0,2) * d_rj(inod,is_fld+2)
!
d_rj(inod,is_rot) = d_rj(inod,is_fld+2)
d_rj(inod,is_rot+1) = d1t_dr1
d_rj(inod,is_rot+2) = - ( d2s_dr2 - g_sph_rj(j,3) &
- & * ar_1d_rj(nlayer_CMB,2)*d_rj(inod,is_fld ) )
+ & * ar_1d_rj(kr_out,2)*d_rj(inod,is_fld ) )
end do
!$omp end parallel do
!
@@ -163,37 +180,38 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_ins_diffuse2(coef_d, &
+ subroutine cal_sph_nod_cmb_ins_diffuse2(jmax, kr_out, &
+ & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, coef_d, &
& is_fld, is_diffuse)
!
- use m_coef_fdm_fixed_CMB
-!
- integer(kind = kint), intent(in) :: is_fld
- integer(kind = kint), intent(in) :: is_diffuse
+ integer(kind = kint), intent(in) :: jmax, kr_out
+ integer(kind = kint), intent(in) :: is_fld, is_diffuse
real(kind = kreal), intent(in) :: coef_d
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_CMB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
!
real(kind = kreal) :: d2s_dr2, d2t_dr2
integer(kind = kint) :: j, inod, i_n1, i_n2
!
!
!$omp parallel do private(inod,i_n1,i_n2,d2s_dr2,d2t_dr2)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_CMB-1) * nidx_rj(2)
- i_n1 = inod - nidx_rj(2)
- i_n2 = i_n1 - nidx_rj(2)
-!
- d2s_dr2 = coef_fdm_fix_dr_CMB_2(-1,3) * d_rj(i_n1,is_fld ) &
- & + (coef_fdm_fix_dr_CMB_2( 0,3) &
- & - g_sph_rj(j,1)*ar_1d_rj(nlayer_CMB,1) &
- & *coef_fdm_fix_dr_CMB_2( 1,3)) * d_rj(inod,is_fld )
- d2t_dr2 = coef_fdm_fix_CMB_2(2,3) * d_rj(i_n2,is_fld+2) &
- & + coef_fdm_fix_CMB_2(1,3) * d_rj(i_n1,is_fld+2) &
- & + coef_fdm_fix_CMB_2(0,3) * d_rj(inod,is_fld+2)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ i_n1 = inod - jmax
+ i_n2 = i_n1 - jmax
+!
+ d2s_dr2 = fdm2_fix_dr_CMB(-1,3) * d_rj(i_n1,is_fld ) &
+ & + (fdm2_fix_dr_CMB( 0,3) &
+ & - g_sph_rj(j,1)*ar_1d_rj(kr_out,1) &
+ & *fdm2_fix_dr_CMB( 1,3)) * d_rj(inod,is_fld )
+ d2t_dr2 = fdm2_fix_fld_CMB(2,3) * d_rj(i_n2,is_fld+2) &
+ & + fdm2_fix_fld_CMB(1,3) * d_rj(i_n1,is_fld+2) &
+ & + fdm2_fix_fld_CMB(0,3) * d_rj(inod,is_fld+2)
!
d_rj(inod,is_diffuse ) = coef_d * (d2s_dr2 &
- & - g_sph_rj(j,3)*ar_1d_rj(nlayer_CMB,2)*d_rj(inod,is_fld ) )
+ & - g_sph_rj(j,3)*ar_1d_rj(kr_out,2)*d_rj(inod,is_fld ) )
d_rj(inod,is_diffuse+2) = coef_d * (d2t_dr2 &
- & - g_sph_rj(j,3)*ar_1d_rj(nlayer_CMB,2)*d_rj(inod,is_fld+2) )
+ & - g_sph_rj(j,3)*ar_1d_rj(kr_out,2)*d_rj(inod,is_fld+2) )
end do
!$omp end parallel do
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_icb_ins.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_icb_ins.f90
index a6ee12c..68f6a1c 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_icb_ins.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_icb_ins.f90
@@ -3,19 +3,31 @@
!!
!!@author H. Matsui
!!@date Programmed in Jan., 2010
+!!@date Modified in Nov., 2013
!
!>@brief Set insulated magnetic boundary condition for ICB
!!
!!@verbatim
-!! subroutine cal_sph_nod_icb_ins_b_and_j(is_fld, is_rot)
-!! subroutine cal_sph_nod_icb_ins_mag2(is_fld)
+!! subroutine cal_sph_nod_icb_ins_b_and_j(jmax, kr_in, &
+!! & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, is_fld, is_rot)
+!! subroutine cal_sph_nod_icb_ins_mag2(jmax, kr_in, is_fld)
!!
-!! subroutine cal_sph_nod_icb_ins_vp_rot2(is_fld, is_rot)
-!! subroutine cal_sph_nod_icb_ins_rot2(is_fld, is_rot)
-!! subroutine cal_sph_nod_icb_ins_diffuse2(coef_d, &
+!! subroutine cal_sph_nod_icb_ins_vp_rot2(jmax, kr_in, &
+!! & is_fld, is_rot)
+!! subroutine cal_sph_nod_icb_ins_rot2(jmax, kr_in, &
+!! & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, is_fld, is_rot)
+!! subroutine cal_sph_nod_icb_ins_diffuse2(jmax, kr_in, &
+!! & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, coef_d, &
!! & is_fld, is_diffuse)
!!@endverbatim
!!
+!!@n @param jmax Number of modes for spherical harmonics @f$L*(L+2)@f$
+!!@n @param kr_in Radial ID for inner boundary
+!!@n @param fdm2_fix_fld_ICB(0:2,3)
+!!! Matrix to evaluate radial derivative at ICB with fiexed field
+!!@n @param fdm2_fix_dr_ICB(-1:1,3)
+!! Matrix to evaluate field at ICB with fiexed radial derivative
+!!
!!@n @param coef_d Coefficient for diffusion term
!!@n @param is_fld Field address of input field
!!@n @param is_rot Field address for curl of field
@@ -40,38 +52,41 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_ins_b_and_j(is_fld, is_rot)
+ subroutine cal_sph_nod_icb_ins_b_and_j(jmax, kr_in, &
+ & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, is_fld, is_rot)
!
use m_coef_fdm_fixed_ICB
!
- integer(kind = kint), intent(in) :: is_fld
- integer(kind = kint), intent(in) :: is_rot
+ integer(kind = kint), intent(in) :: jmax, kr_in
+ integer(kind = kint), intent(in) :: is_fld, is_rot
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_ICB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
!
real(kind = kreal) :: d1s_dr1, d2s_dr2, d1t_dr1
integer(kind = kint) :: j, inod, i_p1, i_p2
!
!
!$omp parallel do private(inod,i_p1,i_p2,d1s_dr1,d2s_dr2,d1t_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_ICB-1) * nidx_rj(2)
- i_p1 = inod + nidx_rj(2)
- i_p2 = i_p1 + nidx_rj(2)
+ do j = 1, jmax
+ inod = j + (kr_in-1) * jmax
+ i_p1 = inod + jmax
+ i_p2 = i_p1 + jmax
!
- d1s_dr1 = (g_sph_rj(j,1)+one) * ar_1d_rj(nlayer_ICB,1) &
+ d1s_dr1 = (g_sph_rj(j,1)+one) * ar_1d_rj(kr_in,1) &
& * d_rj(inod,is_fld)
- d2s_dr2 = ( (g_sph_rj(j,1)+one) * ar_1d_rj(nlayer_ICB,1) &
- & * coef_fdm_fix_dr_ICB_2(-1,3) &
- & + coef_fdm_fix_dr_ICB_2( 0,3) ) * d_rj(inod,is_fld ) &
- & + coef_fdm_fix_dr_ICB_2( 1,3) * d_rj(i_p1,is_fld )
- d1t_dr1 = coef_fdm_fix_ICB_2( 1,2) * d_rj(i_p1,is_fld+2) &
- & + coef_fdm_fix_ICB_2( 2,2) * d_rj(i_p2,is_fld+2)
+ d2s_dr2 = ( (g_sph_rj(j,1)+one) * ar_1d_rj(kr_in,1) &
+ & * fdm2_fix_dr_ICB(-1,3) &
+ & + fdm2_fix_dr_ICB( 0,3) ) * d_rj(inod,is_fld ) &
+ & + fdm2_fix_dr_ICB( 1,3) * d_rj(i_p1,is_fld )
+ d1t_dr1 = fdm2_fix_fld_ICB( 1,2) * d_rj(i_p1,is_fld+2) &
+ & + fdm2_fix_fld_ICB( 2,2) * d_rj(i_p2,is_fld+2)
!
d_rj(inod,is_fld+1) = d1s_dr1
d_rj(inod,is_fld+2) = zero
d_rj(inod,is_rot ) = zero
d_rj(inod,is_rot+1) = d1t_dr1
d_rj(inod,is_rot+2) = - ( d2s_dr2 - g_sph_rj(j,3) &
- & *ar_1d_rj(nlayer_ICB,2)*d_rj(inod,is_fld ) )
+ & *ar_1d_rj(kr_in,2)*d_rj(inod,is_fld ) )
end do
!$omp end parallel do
!
@@ -79,8 +94,9 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_ins_mag2(is_fld)
+ subroutine cal_sph_nod_icb_ins_mag2(jmax, kr_in, is_fld)
!
+ integer(kind = kint), intent(in) :: jmax, kr_in
integer(kind = kint), intent(in) :: is_fld
!
real(kind = kreal) :: d1s_dr1
@@ -88,9 +104,9 @@
!
!
!$omp parallel do private(inod,d1s_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_ICB-1) * nidx_rj(2)
- d1s_dr1 = (g_sph_rj(j,1)+one) * ar_1d_rj(nlayer_ICB,1) &
+ do j = 1, jmax
+ inod = j + (kr_in-1) * jmax
+ d1s_dr1 = (g_sph_rj(j,1)+one) * ar_1d_rj(kr_in,1) &
& * d_rj(inod,is_fld)
!
d_rj(inod,is_fld+1) = d1s_dr1
@@ -103,8 +119,10 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_ins_vp_rot2(is_fld, is_rot)
+ subroutine cal_sph_nod_icb_ins_vp_rot2(jmax, kr_in, &
+ & is_fld, is_rot)
!
+ integer(kind = kint), intent(in) :: jmax, kr_in
integer(kind = kint), intent(in) :: is_fld, is_rot
!
real(kind = kreal) :: d1t_dr1
@@ -112,9 +130,9 @@
!
!
!$omp parallel do private(inod,d1t_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_ICB-1) * nidx_rj(2)
- d1t_dr1 = (g_sph_rj(j,1)+one) * ar_1d_rj(nlayer_ICB,1) &
+ do j = 1, jmax
+ inod = j + (kr_in-1) * jmax
+ d1t_dr1 = (g_sph_rj(j,1)+one) * ar_1d_rj(kr_in,1) &
& * d_rj(inod,is_fld+2)
!
d_rj(inod,is_rot ) = d_rj(inod,is_fld+2)
@@ -127,35 +145,38 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_ins_rot2(is_fld, is_rot)
+ subroutine cal_sph_nod_icb_ins_rot2(jmax, kr_in, &
+ & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, is_fld, is_rot)
!
use m_coef_fdm_fixed_ICB
!
- integer(kind = kint), intent(in) :: is_fld
- integer(kind = kint), intent(in) :: is_rot
+ integer(kind = kint), intent(in) :: jmax, kr_in
+ integer(kind = kint), intent(in) :: is_fld, is_rot
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_ICB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
!
real(kind = kreal) :: d2s_dr2, d1t_dr1
integer(kind = kint) :: j, inod, i_p1, i_p2
!
!
!$omp parallel do private(inod,i_p1,i_p2,d2s_dr2,d1t_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_ICB-1) * nidx_rj(2)
- i_p1 = inod + nidx_rj(2)
- i_p2 = i_p1 + nidx_rj(2)
-!
- d2s_dr2 = ( (g_sph_rj(j,1)+one) * ar_1d_rj(nlayer_ICB,1) &
- & * coef_fdm_fix_dr_ICB_2(-1,3) &
- & + coef_fdm_fix_dr_ICB_2( 0,3) ) * d_rj(inod,is_fld ) &
- & + coef_fdm_fix_dr_ICB_2( 1,3) * d_rj(i_p1,is_fld )
- d1t_dr1 = coef_fdm_fix_ICB_2( 0,2) * d_rj(inod,is_fld+2) &
- & + coef_fdm_fix_ICB_2( 1,2) * d_rj(i_p1,is_fld+2) &
- & + coef_fdm_fix_ICB_2( 2,2) * d_rj(i_p2,is_fld+2)
+ do j = 1, jmax
+ inod = j + (kr_in-1) * jmax
+ i_p1 = inod + jmax
+ i_p2 = i_p1 + jmax
+!
+ d2s_dr2 = ( (g_sph_rj(j,1)+one) * ar_1d_rj(kr_in,1) &
+ & * fdm2_fix_dr_ICB(-1,3) &
+ & + fdm2_fix_dr_ICB( 0,3) ) * d_rj(inod,is_fld ) &
+ & + fdm2_fix_dr_ICB( 1,3) * d_rj(i_p1,is_fld )
+ d1t_dr1 = fdm2_fix_fld_ICB( 0,2) * d_rj(inod,is_fld+2) &
+ & + fdm2_fix_fld_ICB( 1,2) * d_rj(i_p1,is_fld+2) &
+ & + fdm2_fix_fld_ICB( 2,2) * d_rj(i_p2,is_fld+2)
!
d_rj(inod,is_rot ) = d_rj(inod,is_fld+2)
d_rj(inod,is_rot+1) = d1t_dr1
d_rj(inod,is_rot+2) = - ( d2s_dr2 &
- & - g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2)*d_rj(inod,is_fld ) )
+ & - g_sph_rj(j,3)*ar_1d_rj(kr_in,2)*d_rj(inod,is_fld ) )
end do
!$omp end parallel do
!
@@ -163,36 +184,39 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_ins_diffuse2(coef_d, &
+ subroutine cal_sph_nod_icb_ins_diffuse2(jmax, kr_in, &
+ & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, coef_d, &
& is_fld, is_diffuse)
!
use m_coef_fdm_fixed_ICB
!
- integer(kind = kint), intent(in) :: is_fld
- integer(kind = kint), intent(in) :: is_diffuse
+ integer(kind = kint), intent(in) :: jmax, kr_in
+ integer(kind = kint), intent(in) :: is_fld, is_diffuse
real(kind = kreal), intent(in) :: coef_d
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_ICB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
!
real(kind = kreal) :: d2s_dr2,d2t_dr2
integer(kind = kint) :: j, inod,i_p1,i_p2
!
!$omp parallel do private(inod,i_p1,i_p2,d2s_dr2,d2t_dr2)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_ICB-1) * nidx_rj(2)
- i_p1 = inod + nidx_rj(2)
- i_p2 = i_p1 + nidx_rj(2)
-!
- d2s_dr2 = ( (g_sph_rj(j,1)+one) * ar_1d_rj(nlayer_ICB,1) &
- & * coef_fdm_fix_dr_ICB_2(-1,3) &
- & + coef_fdm_fix_dr_ICB_2( 0,3) ) * d_rj(inod,is_fld ) &
- & + coef_fdm_fix_dr_ICB_2( 1,3) * d_rj(i_p1,is_fld )
- d2t_dr2 = coef_fdm_fix_ICB_2( 0,3) * d_rj(inod,is_fld+2) &
- & + coef_fdm_fix_ICB_2( 1,3) * d_rj(i_p1,is_fld+2) &
- & + coef_fdm_fix_ICB_2( 2,3) * d_rj(i_p2,is_fld+2)
+ do j = 1, jmax
+ inod = j + (kr_in-1) * jmax
+ i_p1 = inod + jmax
+ i_p2 = i_p1 + jmax
+!
+ d2s_dr2 = ( (g_sph_rj(j,1)+one) * ar_1d_rj(kr_in,1) &
+ & * fdm2_fix_dr_ICB(-1,3) &
+ & + fdm2_fix_dr_ICB( 0,3) ) * d_rj(inod,is_fld ) &
+ & + fdm2_fix_dr_ICB( 1,3) * d_rj(i_p1,is_fld )
+ d2t_dr2 = fdm2_fix_fld_ICB( 0,3) * d_rj(inod,is_fld+2) &
+ & + fdm2_fix_fld_ICB( 1,3) * d_rj(i_p1,is_fld+2) &
+ & + fdm2_fix_fld_ICB( 2,3) * d_rj(i_p2,is_fld+2)
!
d_rj(inod,is_diffuse ) = coef_d * (d2s_dr2 &
- & - g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2)*d_rj(inod,is_fld ) )
+ & - g_sph_rj(j,3)*ar_1d_rj(kr_in,2)*d_rj(inod,is_fld ) )
d_rj(inod,is_diffuse+2) = coef_d * (d2t_dr2 &
- & - g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2)*d_rj(inod,is_fld+2) )
+ & - g_sph_rj(j,3)*ar_1d_rj(kr_in,2)*d_rj(inod,is_fld+2) )
end do
!$omp end parallel do
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_diffusion.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_diffusion.f90
index 07b8f6d..31dd710 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_diffusion.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_diffusion.f90
@@ -166,8 +166,9 @@
call cal_dsdr_sph_icb_nobc_2(ipol%i_b_diffuse, &
& idpdr%i_b_diffuse)
else
- call cal_sph_nod_icb_ins_diffuse2(coef_d_magne, &
- & ipol%i_magne, ipol%i_b_diffuse)
+ call cal_sph_nod_icb_ins_diffuse2(nidx_rj(2), sph_bc_B%kr_in, &
+ & sph_bc_B%fdm2_fix_fld_ICB, sph_bc_B%fdm2_fix_dr_ICB, &
+ & coef_d_magne, ipol%i_magne, ipol%i_b_diffuse)
call cal_dsdr_sph_icb_nobc_2(ipol%i_b_diffuse, &
& idpdr%i_b_diffuse)
end if
@@ -182,8 +183,9 @@
call cal_sph_nod_cmb_qvc_diffuse2(coef_d_magne, &
& ipol%i_magne, ipol%i_b_diffuse)
else
- call cal_sph_nod_cmb_ins_diffuse2(coef_d_magne, &
- & ipol%i_magne, ipol%i_b_diffuse)
+ call cal_sph_nod_cmb_ins_diffuse2(nidx_rj(2), sph_bc_B%kr_out, &
+ & sph_bc_B%fdm2_fix_fld_CMB, sph_bc_B%fdm2_fix_dr_CMB, &
+ & coef_d_magne, ipol%i_magne, ipol%i_b_diffuse)
end if
call cal_dsdr_sph_cmb_nobc_2(ipol%i_b_diffuse, idpdr%i_b_diffuse)
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_radial_grad.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_radial_grad.f90
index f40aaf4..75cc287 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_radial_grad.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_radial_grad.f90
@@ -170,13 +170,17 @@
call cal_sph_nod_icb_qvc_b_and_j(ipol%i_magne, ipol%i_current)
else
kr_in = nlayer_ICB
- call cal_sph_nod_icb_ins_b_and_j(ipol%i_magne, ipol%i_current)
+ call cal_sph_nod_icb_ins_b_and_j(nidx_rj(2), sph_bc_B%kr_in, &
+ & sph_bc_B%fdm2_fix_fld_ICB, sph_bc_B%fdm2_fix_dr_ICB, &
+ & ipol%i_magne, ipol%i_current)
end if
!
if(sph_bc_B%iflag_cmb .eq. iflag_radial_magne) then
call cal_sph_nod_cmb_qvc_b_and_j(ipol%i_magne, ipol%i_current)
else
- call cal_sph_nod_cmb_ins_b_and_j(ipol%i_magne, ipol%i_current)
+ call cal_sph_nod_cmb_ins_b_and_j(nidx_rj(2), sph_bc_B%kr_out, &
+ & sph_bc_B%fdm2_fix_fld_CMB, sph_bc_B%fdm2_fix_dr_CMB, &
+ & ipol%i_magne, ipol%i_current)
end if
!
!
@@ -259,13 +263,15 @@
call cal_sph_nod_icb_qvc_mag2(ipol%i_magne)
else
kr_in = nlayer_ICB
- call cal_sph_nod_icb_ins_mag2(ipol%i_magne)
+ call cal_sph_nod_icb_ins_mag2(nidx_rj(2), sph_bc_B%kr_in, &
+ & ipol%i_magne)
end if
!
if(sph_bc_B%iflag_cmb .eq. iflag_radial_magne) then
call cal_sph_nod_cmb_qvc_mag2(ipol%i_magne)
else
- call cal_sph_nod_cmb_ins_mag2(ipol%i_magne)
+ call cal_sph_nod_cmb_ins_mag2(nidx_rj(2), sph_bc_B%kr_out, &
+ & ipol%i_magne)
end if
!
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_rotation.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_rotation.f90
index 8382f4c..5a7e748 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_rotation.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_rotation.f90
@@ -116,7 +116,9 @@
call cal_sph_nod_icb_qvc_rot2(ipol%i_magne, ipol%i_current)
else
kr_in = nlayer_ICB
- call cal_sph_nod_icb_ins_rot2(ipol%i_magne, ipol%i_current)
+ call cal_sph_nod_icb_ins_rot2(nidx_rj(2), sph_bc_B%kr_in, &
+ & sph_bc_B%fdm2_fix_fld_ICB, sph_bc_B%fdm2_fix_dr_ICB, &
+ & ipol%i_magne, ipol%i_current)
end if
!
call cal_sph_nod_vect_rot2(kr_in, nlayer_CMB, &
@@ -125,7 +127,9 @@
if(sph_bc_B%iflag_cmb .eq. iflag_radial_magne) then
call cal_sph_nod_cmb_qvc_rot2(ipol%i_magne, ipol%i_current)
else
- call cal_sph_nod_cmb_ins_rot2(ipol%i_magne, ipol%i_current)
+ call cal_sph_nod_cmb_ins_rot2(nidx_rj(2), sph_bc_B%kr_out, &
+ & sph_bc_B%fdm2_fix_fld_CMB, sph_bc_B%fdm2_fix_dr_CMB, &
+ & ipol%i_magne, ipol%i_current)
end if
!
end subroutine const_sph_current
@@ -155,7 +159,8 @@
call cal_sph_nod_icb_qvc_vp_rot2(is_fld, is_rot)
else
kr_st = nlayer_ICB
- call cal_sph_nod_icb_ins_vp_rot2(is_fld, is_rot)
+ call cal_sph_nod_icb_ins_vp_rot2(nidx_rj(2), sph_bc_B%kr_in, &
+ & is_fld, is_rot)
end if
!
call cal_sph_nod_vect_w_div_rot2(kr_st, nlayer_CMB, &
@@ -164,7 +169,8 @@
if(sph_bc_B%iflag_cmb .eq. iflag_radial_magne) then
call cal_sph_nod_cmb_qvc_vp_rot2(is_fld, is_rot)
else
- call cal_sph_nod_cmb_ins_vp_rot2(is_fld, is_rot)
+ call cal_sph_nod_cmb_ins_vp_rot2(nidx_rj(2), sph_bc_B%kr_out, &
+ & is_fld, is_rot)
end if
!
end subroutine const_sph_rotation_uxb
@@ -297,8 +303,9 @@
& idpdr%i_b_diffuse)
else
kr_in = nlayer_ICB
- call cal_sph_nod_icb_ins_diffuse2(coef_d_magne, &
- & ipol%i_magne, ipol%i_b_diffuse)
+ call cal_sph_nod_icb_ins_diffuse2(nidx_rj(2), sph_bc_B%kr_in, &
+ & sph_bc_B%fdm2_fix_fld_ICB, sph_bc_B%fdm2_fix_dr_ICB, &
+ & coef_d_magne, ipol%i_magne, ipol%i_b_diffuse)
call cal_dsdr_sph_icb_nobc_2(ipol%i_b_diffuse, &
& idpdr%i_b_diffuse)
end if
@@ -310,8 +317,9 @@
call cal_sph_nod_cmb_qvc_diffuse2(coef_d_magne, &
& ipol%i_magne, ipol%i_b_diffuse)
else
- call cal_sph_nod_cmb_ins_diffuse2(coef_d_magne, &
- & ipol%i_magne, ipol%i_b_diffuse)
+ call cal_sph_nod_cmb_ins_diffuse2(nidx_rj(2), sph_bc_B%kr_out, &
+ & sph_bc_B%fdm2_fix_fld_CMB, sph_bc_B%fdm2_fix_dr_CMB, &
+ & coef_d_magne, ipol%i_magne, ipol%i_b_diffuse)
end if
call cal_dsdr_sph_cmb_nobc_2(ipol%i_b_diffuse, idpdr%i_b_diffuse)
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/m_addresses_trans_sph_MHD.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/m_addresses_trans_sph_MHD.f90
index 79bbf0c..8635f21 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/m_addresses_trans_sph_MHD.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/m_addresses_trans_sph_MHD.f90
@@ -72,13 +72,13 @@
b_trns%i_current = nvector_rj_2_rtp
end if
!
-! gradient of temperature flag
+! temperature flag
if(iflag_t_evo_4_temp .gt. id_no_evolution) then
nscalar_rj_2_rtp = nscalar_rj_2_rtp + 1
b_trns%i_temp = nscalar_rj_2_rtp
end if
!
-! gradient of dummy scalar flag
+! composition flag
if(iflag_t_evo_4_composit .gt. id_no_evolution) then
nscalar_rj_2_rtp = nscalar_rj_2_rtp + 1
b_trns%i_light = nscalar_rj_2_rtp
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_bc_sph_mhd.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_bc_sph_mhd.f90
index 4c422a2..c71c708 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_bc_sph_mhd.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_bc_sph_mhd.f90
@@ -37,7 +37,13 @@
use set_bc_flag_sph_velo
use set_bc_sph_scalars
use set_reference_sph_mhd
- use cal_sph_bc_fdm_matrix
+!
+ use m_coef_fdm_fixed_ICB
+ use m_coef_fdm_fixed_CMB
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
+ use m_coef_fdm_to_center
+ use cal_fdm_coefs_4_boundaries
!
!
if (iflag_t_evo_4_velo .gt. id_no_evolution) then
@@ -57,10 +63,40 @@
end if
!
!
+! Set FDM matrices for boundaries
+!
call set_radial_range_by_BC(sph_bc_U)
+ call set_radial_range_by_BC(sph_bc_B)
call set_radial_range_by_BC(sph_bc_T)
call set_radial_range_by_BC(sph_bc_C)
- call set_radial_range_by_BC(sph_bc_B)
+!
+! Set FDM matrices for boundaries
+!
+ call cal_fdm_coefs_4_BCs(nidx_rj(1), radius_1d_rj_r, sph_bc_U)
+ call cal_fdm_coefs_4_BCs(nidx_rj(1), radius_1d_rj_r, sph_bc_B)
+ call cal_fdm_coefs_4_BCs(nidx_rj(1), radius_1d_rj_r, sph_bc_T)
+ call cal_fdm_coefs_4_BCs(nidx_rj(1), radius_1d_rj_r, sph_bc_C)
+!
+ call cal_fdm2_coef_fix_fld_ICB(radius_1d_rj_r(nlayer_ICB), &
+ & coef_fdm_fix_ICB_2)
+ call cal_fdm2_coef_fix_df_ICB(radius_1d_rj_r(nlayer_ICB), &
+ & coef_fdm_fix_dr_ICB_2)
+!
+ call cal_fdm2_coef_fix_fld_CMB(radius_1d_rj_r(nlayer_CMB-2), &
+ & coef_fdm_fix_CMB_2)
+ call cal_fdm2_coef_fix_df_CMB(radius_1d_rj_r(nlayer_CMB-1), &
+ & coef_fdm_fix_dr_CMB_2)
+!
+!
+ call cal_2nd_ICB_free_vp_bc_fdm(radius_1d_rj_r(nlayer_ICB))
+ call cal_2nd_ICB_free_vt_bc_fdm(radius_1d_rj_r(nlayer_ICB))
+!
+ call cal_2nd_CMB_free_vp_bc_fdm(radius_1d_rj_r(nlayer_CMB-1))
+ call cal_2nd_CMB_free_vt_bc_fdm(radius_1d_rj_r(nlayer_CMB-1))
+!
+ call cal_2nd_to_center_fixed_fdm(radius_1d_rj_r(1))
+ call cal_2nd_to_center_fix_df_fdm(radius_1d_rj_r(1))
+!
!
! Set reference temperature and adjust boundary conditions
!
@@ -68,20 +104,28 @@
call set_ref_temp_sph_mhd
call adjust_sph_temp_bc_by_reftemp
!
-! Det FDM matrices for boundaries
-!
- if (iflag_debug.gt.0) write(*,*) 's_cal_sph_bc_fdm_matrices'
- call s_cal_sph_bc_fdm_matrices
-!
! Check data
!
if(i_debug .gt. 1) then
- call check_sph_boundary_spectra(fhd_temp, &
- & nidx_rj(2), idx_gl_1d_rj_j, sph_bc_T)
+ if (iflag_t_evo_4_temp .gt. id_no_evolution) then
+ call check_sph_boundary_spectra(fhd_temp, &
+ & nidx_rj(2), idx_gl_1d_rj_j, sph_bc_T)
+ end if
+ if (iflag_t_evo_4_composit .gt. id_no_evolution) then
+ call check_sph_boundary_spectra(fhd_light, &
+ & nidx_rj(2), idx_gl_1d_rj_j, sph_bc_C)
+ end if
end if
- if(i_debug .gt. 1) then
- call check_sph_boundary_spectra(fhd_light, &
- & nidx_rj(2), idx_gl_1d_rj_j, sph_bc_C)
+!
+ if (iflag_debug .eq. iflag_full_msg) then
+ call check_fdm_coefs_4_BC2(fhd_velo, sph_bc_U)
+ call check_fdm_coefs_4_BC2(fhd_magne, sph_bc_B)
+ call check_fdm_coefs_4_BC2(fhd_temp, sph_bc_T)
+ call check_fdm_coefs_4_BC2(fhd_light, sph_bc_C)
+ call check_coef_fdm_fix_dr_CMB
+ call check_coef_fdm_free_ICB
+ call check_coef_fdm_free_CMB
+ call check_coef_fdm_fix_dr_2ctr
end if
!
end subroutine s_set_bc_sph_mhd
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90
index 15e3f5f..4fb91c3 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90
@@ -122,9 +122,9 @@
& sph_bc%fdm2_fix_dr_ICB)
!
call cal_fdm2_coef_fix_fld_CMB(radius(sph_bc%kr_out-2), &
- & sph_bc%fdm2_fix_fld_ICB)
+ & sph_bc%fdm2_fix_fld_CMB)
call cal_fdm2_coef_fix_df_CMB(radius(sph_bc%kr_out-1), &
- & sph_bc%fdm2_fix_dr_ICB)
+ & sph_bc%fdm2_fix_dr_CMB)
!
end subroutine cal_fdm_coefs_4_BCs
!
More information about the CIG-COMMITS
mailing list