[cig-commits] [commit] Hiro_latest: Use structures for free slip boundaries (59a80c8)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Nov 18 16:22:03 PST 2013
Repository : ssh://geoshell/calypso
On branch : Hiro_latest
Link : https://github.com/geodynamics/calypso/compare/93e9f8f974c7a247c8f02e54ec18de063f86c8fb...3c548304673360ddedd7d68c8095b3fb74a2b9ce
>---------------------------------------------------------------
commit 59a80c8350553f0d91dd0351d46d20c8624c314b
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Fri Nov 15 11:28:33 2013 -0800
Use structures for free slip boundaries
>---------------------------------------------------------------
59a80c8350553f0d91dd0351d46d20c8624c314b
.../MHD_src/sph_MHD/cal_sol_sph_fluid_crank.f90 | 8 +-
.../MHD_src/sph_MHD/cal_sph_exp_fixed_flux.f90 | 20 ++-
.../MHD_src/sph_MHD/const_sph_diffusion.f90 | 28 ++--
.../MHD_src/sph_MHD/const_sph_radial_grad.f90 | 18 ++-
.../MHD_src/sph_MHD/const_sph_rotation.f90 | 38 ++++--
.../MHD_src/sph_MHD/set_sph_exp_free_CMB.f90 | 148 +++++++++++++--------
.../MHD_src/sph_MHD/set_sph_exp_free_ICB.f90 | 148 +++++++++++++--------
.../MHD_src/sph_MHD/set_sph_exp_rigid_ICB.f90 | 12 ++
8 files changed, 278 insertions(+), 142 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 e5af4f7..4f3afec 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
@@ -53,6 +53,8 @@
subroutine cal_sol_velo_by_vort_sph_crank
!
use m_boundary_params_sph_MHD
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
use const_sph_radial_grad
use set_sph_exp_rigid_ICB
use set_sph_exp_rigid_CMB
@@ -73,7 +75,8 @@
call delete_zero_degree_comp(ipol%i_velo)
!
if (sph_bc_U%iflag_icb .eq. iflag_free_slip) then
- call cal_sph_nod_icb_free_vpol2(ipol%i_velo)
+ call cal_sph_nod_icb_free_vpol2(nidx_rj(2), sph_bc_U%kr_in, &
+ & coef_fdm_free_ICB_vp2, ipol%i_velo)
else if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then
call cal_sph_nod_icb_rotate_velo2 &
& (idx_rj_degree_zero, idx_rj_degree_one, nidx_rj(2), &
@@ -84,7 +87,8 @@
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
- call cal_sph_nod_cmb_free_vpol2(ipol%i_velo)
+ call cal_sph_nod_cmb_free_vpol2(nidx_rj(2), sph_bc_U%kr_out, &
+ & coef_fdm_free_CMB_vp2, ipol%i_velo)
else
call cal_sph_nod_cmb_rigid_velo2(nidx_rj(2), &
& sph_bc_U%kr_out, sph_bc_U%r_CMB, vt_CMB_bc, ipol%i_velo)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_fixed_flux.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_fixed_flux.f90
index 89bf4a5..4e5d6b8 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_fixed_flux.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_fixed_flux.f90
@@ -12,7 +12,7 @@
!! subroutine cal_div_sph_in_fix_flux_2(jmax, kr_in, r_ICB, &
!! & flux_ICB, is_fld, is_div)
!! subroutine cal_sph_in_fix_flux_diffuse2(jmax, kr_in, r_ICB, &
-!! & fdm2_fix_dr_ICB, flux_IN, coef_d, is_fld, is_diffuse)
+!! & fdm2_fix_dr_ICB, flux_ICB, coef_d, is_fld, is_diffuse)
!!
!! subroutine cal_dsdr_sph_out_fix_flux_2(idx_rj_degree_zero, &
!! & jmax, kr_out, r_CMB, flux_CMB, is_fld, is_grd)
@@ -22,9 +22,19 @@
!! & fdm2_fix_dr_CMB, flux_OUT, coef_d, is_fld, is_diffuse)
!!@endverbatim
!!
+!!@n @param idx_rj_degree_zero Local address for degree 0
!!@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 kr_out Radial ID for outer boundary
+!!@n @param r_ICB(0:2) Radius at ICB
+!!@n @param r_CMB(0:2) Radius at CMB
!!@n @param flux_ICB(jamx) Spectrum of fixed flux at ICB
!!@n @param flux_CMB(jamx) Spectrum of fixed flux at CMB
+!!@n @param fdm2_fix_dr_ICB(-1:1,3)
+!! Matrix to evaluate field at ICB with fiexed radial derivative
+!!@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 Address of spectrum data d_rj
@@ -104,12 +114,12 @@
! -----------------------------------------------------------------------
!
subroutine cal_sph_in_fix_flux_diffuse2(jmax, kr_in, r_ICB, &
- & fdm2_fix_dr_ICB, flux_IN, coef_d, is_fld, is_diffuse)
+ & fdm2_fix_dr_ICB, flux_ICB, coef_d, is_fld, 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) :: flux_IN(jmax)
+ real(kind = kreal), intent(in) :: flux_ICB(jmax)
real(kind = kreal), intent(in) :: r_ICB(0:2)
real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
!
@@ -122,12 +132,12 @@
inod = j + (kr_in-1) * jmax
i_p1 = inod + jmax
!
- d2t_dr2 = fdm2_fix_dr_ICB(-1,3) * flux_IN(j) &
+ d2t_dr2 = fdm2_fix_dr_ICB(-1,3) * flux_ICB(j) &
& + fdm2_fix_dr_ICB( 0,3) * d_rj(inod,is_fld) &
& + fdm2_fix_dr_ICB( 1,3) * d_rj(i_p1,is_fld)
!
d_rj(inod,is_diffuse) = coef_d * (d2t_dr2 &
- & + two*r_ICB(1) * flux_IN(j) &
+ & + two*r_ICB(1) * flux_ICB(j) &
& - g_sph_rj(j,3)*r_ICB(2) * d_rj(inod,is_fld))
!
end 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 d88c010..400f37f 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
@@ -47,6 +47,8 @@
subroutine const_sph_viscous_diffusion
!
use m_boundary_params_sph_MHD
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
use set_sph_exp_rigid_ICB
use set_sph_exp_rigid_CMB
use set_sph_exp_free_ICB
@@ -66,8 +68,10 @@
& ipol%i_velo, ipol%i_v_diffuse)
call cal_dsdr_sph_center_2(ipol%i_v_diffuse)
else if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then
- call cal_sph_nod_icb_free_diffuse2(coef_d_velo, &
- & ipol%i_velo, ipol%i_v_diffuse)
+ call cal_sph_nod_icb_free_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+ & coef_d_velo, ipol%i_velo, ipol%i_v_diffuse)
else
call cal_sph_nod_icb_rigid_diffuse2 &
& (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
@@ -85,8 +89,10 @@
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
- call cal_sph_nod_cmb_free_diffuse2(coef_d_velo, &
- & ipol%i_velo, ipol%i_v_diffuse)
+ call cal_sph_nod_cmb_free_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+ & coef_d_velo, ipol%i_velo, ipol%i_v_diffuse)
else
call cal_sph_nod_cmb_rigid_diffuse2 &
& (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
@@ -104,6 +110,8 @@
subroutine const_sph_vorticirty_diffusion
!
use m_boundary_params_sph_MHD
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
use set_sph_exp_rigid_ICB
use set_sph_exp_rigid_CMB
use set_sph_exp_free_ICB
@@ -124,8 +132,10 @@
& ipol%i_vort, ipol%i_w_diffuse)
call cal_dsdr_sph_center_2(ipol%i_w_diffuse)
else if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then
- call cal_sph_nod_icb_free_w_diffuse2(coef_d_velo, &
- & ipol%i_vort, ipol%i_w_diffuse)
+ call cal_sph_nod_icb_free_w_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & sph_bc_U%fdm2_fix_fld_ICB, coef_fdm_free_ICB_vt2, &
+ & coef_d_velo, ipol%i_vort, ipol%i_w_diffuse)
else
call cal_sph_nod_icb_rgd_w_diffuse2(nidx_rj(2), &
& sph_bc_U%kr_in, sph_bc_U%r_ICB, sph_bc_U%fdm2_fix_fld_ICB, &
@@ -142,8 +152,10 @@
& ipol%i_w_diffuse, idpdr%i_w_diffuse)
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
- call cal_sph_nod_cmb_free_w_diffuse2(coef_d_velo, &
- & ipol%i_vort, ipol%i_w_diffuse)
+ call cal_sph_nod_cmb_free_w_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & sph_bc_U%fdm2_fix_fld_CMB, coef_fdm_free_CMB_vt2, &
+ & coef_d_velo, ipol%i_vort, ipol%i_w_diffuse)
else
call cal_sph_nod_cmb_rgd_w_diffuse2(nidx_rj(2), &
& sph_bc_U%kr_out, sph_bc_U%r_CMB, sph_bc_U%fdm2_fix_fld_CMB, &
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 605a871..4e77244 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
@@ -130,6 +130,8 @@
subroutine const_grad_vp_and_vorticity
!
use m_boundary_params_sph_MHD
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
use set_sph_exp_rigid_ICB
use set_sph_exp_rigid_CMB
use set_sph_exp_free_ICB
@@ -138,7 +140,9 @@
!
!
if (sph_bc_U%iflag_icb .eq. iflag_free_slip) then
- call cal_sph_nod_icb_free_v_and_w(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_icb_free_v_and_w(nidx_rj(2), sph_bc_U%kr_in, &
+ & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+ & ipol%i_velo, ipol%i_vort)
else if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then
call cal_sph_nod_icb_rotate_velo2 &
& (idx_rj_degree_zero, idx_rj_degree_one, nidx_rj(2), &
@@ -157,7 +161,9 @@
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
- call cal_sph_nod_cmb_free_v_and_w(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_cmb_free_v_and_w(nidx_rj(2), sph_bc_U%kr_out, &
+ & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+ & ipol%i_velo, ipol%i_vort)
else
call cal_sph_nod_cmb_rigid_v_and_w &
& (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
@@ -235,6 +241,8 @@
subroutine const_grad_poloidal_moment(i_field)
!
use m_boundary_params_sph_MHD
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
use set_sph_exp_rigid_ICB
use set_sph_exp_rigid_CMB
use set_sph_exp_free_ICB
@@ -245,7 +253,8 @@
!
!
if (sph_bc_U%iflag_icb .eq. iflag_free_slip) then
- call cal_sph_nod_icb_free_vpol2(i_field)
+ call cal_sph_nod_icb_free_vpol2(nidx_rj(2), sph_bc_U%kr_in, &
+ & coef_fdm_free_ICB_vp2, i_field)
else if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then
call cal_sph_nod_icb_rotate_velo2 &
& (idx_rj_degree_zero, idx_rj_degree_one, nidx_rj(2), &
@@ -256,7 +265,8 @@
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
- call cal_sph_nod_cmb_free_vpol2(i_field)
+ call cal_sph_nod_cmb_free_vpol2(nidx_rj(2), sph_bc_U%kr_out, &
+ & coef_fdm_free_CMB_vp2, i_field)
else
call cal_sph_nod_cmb_rigid_velo2(nidx_rj(2), &
& sph_bc_U%kr_out, sph_bc_U%r_CMB, vt_CMB_bc, i_field)
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 39f77f0..6a6f94c 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
@@ -67,6 +67,8 @@
subroutine const_sph_vorticity
!
use m_boundary_params_sph_MHD
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
use m_sph_phys_address
use set_sph_exp_rigid_ICB
use set_sph_exp_rigid_CMB
@@ -75,7 +77,10 @@
!
!
if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then
- call cal_sph_nod_icb_free_rot2(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_icb_free_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+ & ipol%i_velo, ipol%i_vort)
else if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then
call cal_sph_nod_icb_rigid_rot2 &
& (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
@@ -89,7 +94,10 @@
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
- call cal_sph_nod_cmb_free_rot2(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_cmb_free_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+ & ipol%i_velo, ipol%i_vort)
else
call cal_sph_nod_cmb_rigid_rot2 &
& (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
@@ -221,6 +229,8 @@
subroutine const_sph_force_rot2(is_fld, is_rot)
!
use m_boundary_params_sph_MHD
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
use set_sph_exp_rigid_ICB
use set_sph_exp_rigid_CMB
use set_sph_exp_free_ICB
@@ -230,7 +240,10 @@
!
!
if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then
- call cal_sph_nod_icb_free_rot2(is_fld, is_rot)
+ call cal_sph_nod_icb_free_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+ & is_fld, is_rot)
else
call cal_sph_nod_icb_rigid_rot2 &
& (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
@@ -239,7 +252,10 @@
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
- call cal_sph_nod_cmb_free_rot2(is_fld, is_rot)
+ call cal_sph_nod_cmb_free_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+ & is_fld, is_rot)
else
call cal_sph_nod_cmb_rigid_rot2 &
& (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
@@ -258,6 +274,8 @@
subroutine const_sph_viscous_by_vort2
!
use m_boundary_params_sph_MHD
+ use m_coef_fdm_free_ICB
+ use m_coef_fdm_free_CMB
use m_sph_phys_address
use m_physical_property
use set_sph_exp_rigid_ICB
@@ -272,8 +290,10 @@
& coef_d_velo, ipol%i_vort, ipol%i_v_diffuse)
!
if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then
- call cal_sph_nod_icb_free_diffuse2(coef_d_velo, &
- & ipol%i_velo, ipol%i_v_diffuse)
+ call cal_sph_nod_icb_free_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+ & coef_d_velo, ipol%i_velo, ipol%i_v_diffuse)
else
call cal_sph_nod_icb_rigid_diffuse2 &
& (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
@@ -290,8 +310,10 @@
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
- call cal_sph_nod_cmb_free_diffuse2(coef_d_velo, &
- & ipol%i_velo, ipol%i_v_diffuse)
+ call cal_sph_nod_cmb_free_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+ & coef_d_velo, ipol%i_velo, ipol%i_v_diffuse)
else
call cal_sph_nod_cmb_rigid_diffuse2 &
& (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_free_CMB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_free_CMB.f90
index 6ba43e8..0337ed0 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_free_CMB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_free_CMB.f90
@@ -7,15 +7,35 @@
!>@brief Evaluate velocity with free slip boundary at CMB
!!
!!@verbatim
-!! subroutine cal_sph_nod_cmb_free_v_and_w(is_fld, is_rot)
-!! subroutine cal_sph_nod_cmb_free_vpol2(is_fld)
-!! subroutine cal_sph_nod_cmb_free_rot2(is_fld, is_rot)
-!! subroutine cal_sph_nod_cmb_free_diffuse2(coef_d, &
-!! & is_fld, is_diffuse)
-!! subroutine cal_sph_nod_cmb_free_w_diffuse2(coef_d, &
-!! & is_fld, is_diffuse)
+!! subroutine cal_sph_nod_cmb_free_v_and_w(jmax, kr_out, &
+!! & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+!! & is_fld, is_rot)
+!! subroutine cal_sph_nod_cmb_free_vpol2(jmax, kr_out, &
+!! & coef_fdm_free_CMB_vp2, is_fld)
+!! subroutine cal_sph_nod_cmb_free_rot2(jmax, kr_out, r_CMB, &
+!! & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+!! & is_fld, is_rot)
+!! subroutine cal_sph_nod_cmb_free_diffuse2(jmax, kr_out, r_CMB, &
+!! & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+!! & coef_d, is_fld, is_diffuse)
+!! subroutine cal_sph_nod_cmb_free_w_diffuse2(jmax, kr_out, r_CMB, &
+!! & fdm2_fix_fld_CMB, coef_fdm_free_CMB_vt2, &
+!! & coef_d, is_fld, is_diffuse)
!!@endverbatim
!!
+!!@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 r_CMB(0:2) Radius at CMB
+!!
+!!@n @param fdm2_fix_fld_CMB(0:2,3)
+!! Matrix to evaluate radial derivative at CMB with fiexed field
+!!@n @param coef_fdm_free_CMB_vp2(-1:0,3)
+!! Matrix to evaluate poloidal velocity
+!! with free slip boundary at CMB
+!!@n @param coef_fdm_free_CMB_vt2(-1:0,3)
+!! Matrix to evaluate toroidal velocity
+!! with free slip boundary at CMB
+!!
!!@n @param coef_d Coefficient for diffusion term
!!@n @param is_fld Address of poloidal velocity in d_rj
!!@n @param is_rot Address of poloidal vorticity in d_rj
@@ -26,10 +46,8 @@
use m_precision
!
use m_constants
- use m_spheric_parameter
use m_schmidt_poly_on_rtm
use m_sph_spectr_data
- use m_coef_fdm_free_CMB
!
implicit none
!
@@ -39,20 +57,25 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_free_v_and_w(is_fld, is_rot)
+ subroutine cal_sph_nod_cmb_free_v_and_w(jmax, kr_out, &
+ & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+ & 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) :: coef_fdm_free_CMB_vp2(-1:0,3)
+ real(kind = kreal), intent(in) :: coef_fdm_free_CMB_vt2(-1:0,3)
!
real(kind = kreal) :: d1s_dr1, d2s_dr2, d1t_dr1
integer(kind = kint) :: inod, j, 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 = coef_fdm_free_CMB_vp2(-1,2) * d_rj(i_n1,is_fld )
d2s_dr2 = coef_fdm_free_CMB_vp2(-1,3) * d_rj(i_n1,is_fld )
@@ -71,18 +94,21 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_free_vpol2(is_fld)
+ subroutine cal_sph_nod_cmb_free_vpol2(jmax, kr_out, &
+ & coef_fdm_free_CMB_vp2, is_fld)
!
+ integer(kind = kint), intent(in) :: jmax, kr_out
integer(kind = kint), intent(in) :: is_fld
+ real(kind = kreal), intent(in) :: coef_fdm_free_CMB_vp2(-1:0,3)
!
real(kind = kreal) :: d1s_dr1
integer(kind = kint) :: inod, j, i_n1
!
!
!$omp parallel do private(inod,i_n1,j,d1s_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_CMB-1) * nidx_rj(2)
- i_n1 = inod - nidx_rj(2)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ i_n1 = inod - jmax
!
d1s_dr1 = coef_fdm_free_CMB_vp2(-1,2) * d_rj(i_n1,is_fld ) &
& + coef_fdm_free_CMB_vp2( 0,2) * d_rj(inod,is_fld )
@@ -96,20 +122,26 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_free_rot2(is_fld, is_rot)
+ subroutine cal_sph_nod_cmb_free_rot2(jmax, kr_out, r_CMB, &
+ & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+ & 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) :: r_CMB(0:2)
+!
+ real(kind = kreal), intent(in) :: coef_fdm_free_CMB_vp2(-1:0,3)
+ real(kind = kreal), intent(in) :: coef_fdm_free_CMB_vt2(-1:0,3)
!
real(kind = kreal) :: d2s_dr2, d1t_dr1
integer(kind = kint) :: inod, j, 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)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ i_n1 = inod - jmax
+ i_n2 = i_n1 - jmax
!
d2s_dr2 = coef_fdm_free_CMB_vp2(-1,3) * d_rj(i_n1,is_fld ) &
& + coef_fdm_free_CMB_vp2( 0,3) * d_rj(inod,is_fld )
@@ -119,8 +151,7 @@
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 )
+ & + g_sph_rj(j,3)*r_CMB(2) * d_rj(inod,is_fld )
end do
!$omp end parallel do
!
@@ -128,21 +159,26 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_free_diffuse2(coef_d, &
- & is_fld, is_diffuse)
+ subroutine cal_sph_nod_cmb_free_diffuse2(jmax, kr_out, r_CMB, &
+ & coef_fdm_free_CMB_vp2, coef_fdm_free_CMB_vt2, &
+ & coef_d, is_fld, is_diffuse)
!
- 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) :: r_CMB(0:2)
real(kind = kreal), intent(in) :: coef_d
!
+ real(kind = kreal), intent(in) :: coef_fdm_free_CMB_vp2(-1:0,3)
+ real(kind = kreal), intent(in) :: coef_fdm_free_CMB_vt2(-1:0,3)
+!
real(kind = kreal) :: d2s_dr2, d2t_dr2
integer(kind = kint) :: inod, j, i_n1
!
!
!$omp parallel do private(inod,i_n1,j,d2s_dr2,d2t_dr2)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_CMB-1) * nidx_rj(2)
- i_n1 = inod - nidx_rj(2)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ i_n1 = inod - jmax
!
d2s_dr2 = coef_fdm_free_CMB_vp2(-1,3) * d_rj(i_n1,is_fld ) &
& + coef_fdm_free_CMB_vp2( 0,3) * d_rj(inod,is_fld )
@@ -150,11 +186,9 @@
& + coef_fdm_free_CMB_vt2( 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)*r_CMB(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)*r_CMB(2) * d_rj(inod,is_fld+2) )
end do
!$omp end parallel do
!
@@ -163,37 +197,37 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_free_w_diffuse2(coef_d, &
- & is_fld, is_diffuse)
+ subroutine cal_sph_nod_cmb_free_w_diffuse2(jmax, kr_out, r_CMB, &
+ & fdm2_fix_fld_CMB, coef_fdm_free_CMB_vt2, &
+ & 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) :: r_CMB(0:2)
real(kind = kreal), intent(in) :: coef_d
+ real(kind = kreal), intent(in) :: coef_fdm_free_CMB_vt2(-1:0,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_CMB(0:2,3)
!
integer(kind = kint) :: inod, j, i_n1, i_n2
real(kind = kreal) :: d2s_dr2, d2t_dr2
!
!
!$omp parallel do private(inod,i_n1,i_n2,j,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)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ i_n1 = inod - jmax
+ i_n2 = i_n1 - jmax
!
d2s_dr2 = coef_fdm_free_CMB_vt2(-1,3) * d_rj(i_n1,is_fld ) &
& + coef_fdm_free_CMB_vt2( 0,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)
+ 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)*r_CMB(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)*r_CMB(2) * d_rj(inod,is_fld+2) )
end do
!$omp end parallel do
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_free_ICB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_free_ICB.f90
index 6e7e38e..ee9b389 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_free_ICB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_free_ICB.f90
@@ -7,15 +7,37 @@
!>@brief Evaluate velocity with free slip boundary at ICB
!!
!!@verbatim
-!! subroutine cal_sph_nod_icb_free_v_and_w(is_fld, is_rot)
-!! subroutine cal_sph_nod_icb_free_vpol2(is_fld)
-!! subroutine cal_sph_nod_icb_free_rot2(is_fld, is_rot)
-!! subroutine cal_sph_nod_icb_free_diffuse2(coef_d, &
-!! & is_fld, is_diffuse)
-!! subroutine cal_sph_nod_icb_free_w_diffuse2(coef_d, &
-!! & is_fld, is_diffuse)
+!! subroutine cal_sph_nod_icb_free_v_and_w(jmax, kr_in, &
+!! & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+!! & is_fld, is_rot)
+!! subroutine cal_sph_nod_icb_free_vpol2(jmax, kr_in, &
+!! & coef_fdm_free_ICB_vp2, is_fld)
+!! subroutine cal_sph_nod_icb_free_rot2(jmax, kr_in, r_ICB, &
+!! & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+!! & is_fld, is_rot)
+!! subroutine cal_sph_nod_icb_free_diffuse2(jmax, kr_in, r_ICB, &
+!! & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+!! & coef_d, is_fld, is_diffuse)
+!! subroutine cal_sph_nod_icb_free_w_diffuse2(jmax, kr_in, r_ICB, &
+!! & fdm2_fix_fld_ICB, coef_fdm_free_ICB_vt2, &
+!! & coef_d, is_fld, is_diffuse)
!!@endverbatim
!!
+!!@n @param idx_rj_degree_zero Local address for degree 0
+!!@n @param idx_rj_degree_one(-1:1) Local address for degree 1
+!!@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 r_ICB(0:2) Radius at ICB
+!!
+!!@n @param fdm2_fix_fld_ICB(0:2,3)
+!! Matrix to evaluate radial derivative at ICB with fiexed field
+!!@n @param coef_fdm_free_ICB_vp2(0:1,3)
+!! Matrix to evaluate poloidal velocity
+!! with free slip boundary at ICB
+!!@n @param coef_fdm_free_ICB_vt2(0:1,3)
+!! Matrix to evaluate toroidal velocity
+!! with free slip boundary at ICB
+!!
!!@n @param coef_d Coefficient for diffusion term
!!@n @param is_fld Address of poloidal velocity in d_rj
!!@n @param is_rot Address of poloidal vorticity in d_rj
@@ -26,10 +48,8 @@
use m_precision
!
use m_constants
- use m_spheric_parameter
use m_schmidt_poly_on_rtm
use m_sph_spectr_data
- use m_coef_fdm_free_ICB
!
implicit none
!
@@ -39,21 +59,24 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_free_v_and_w(is_fld, is_rot)
-!
- integer(kind = kint), intent(in) :: is_fld
- integer(kind = kint), intent(in) :: is_rot
+ subroutine cal_sph_nod_icb_free_v_and_w(jmax, kr_in, &
+ & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+ & is_fld, is_rot)
!
+ integer(kind = kint), intent(in) :: jmax, kr_in
+ integer(kind = kint), intent(in) :: is_fld, is_rot
+ real(kind = kreal), intent(in) :: coef_fdm_free_ICB_vp2(0:1,3)
+ real(kind = kreal), intent(in) :: coef_fdm_free_ICB_vt2(0:1,3)
!
real(kind = kreal) :: d1s_dr1, d2s_dr2, d1t_dr1
integer(kind = kint) :: inod, j, i_p1, i_p2
!
!
!$omp parallel do private(inod,i_p1,i_p2,j,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 = coef_fdm_free_ICB_vp2( 1,2) * d_rj(i_p1,is_fld )
d2s_dr2 = coef_fdm_free_ICB_vp2( 1,3) * d_rj(i_p1,is_fld )
@@ -71,9 +94,12 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_free_vpol2(is_fld)
+ subroutine cal_sph_nod_icb_free_vpol2(jmax, kr_in, &
+ & coef_fdm_free_ICB_vp2, is_fld)
!
+ integer(kind = kint), intent(in) :: jmax, kr_in
integer(kind = kint), intent(in) :: is_fld
+ real(kind = kreal), intent(in) :: coef_fdm_free_ICB_vp2(0:1,3)
!
!
real(kind = kreal) :: d1s_dr1
@@ -81,9 +107,9 @@
!
!
!$omp parallel do private(inod,i_p1,j,d1s_dr1)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_ICB-1) * nidx_rj(2)
- i_p1 = inod + nidx_rj(2)
+ do j = 1, jmax
+ inod = j + (kr_in-1) * jmax
+ i_p1 = inod + jmax
!
d1s_dr1 = coef_fdm_free_ICB_vp2( 0,2) * d_rj(inod,is_fld ) &
& + coef_fdm_free_ICB_vp2( 1,2) * d_rj(i_p1,is_fld )
@@ -97,10 +123,15 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_free_rot2(is_fld, is_rot)
+ subroutine cal_sph_nod_icb_free_rot2(jmax, kr_in, r_ICB, &
+ & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+ & 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_in
+ integer(kind = kint), intent(in) :: is_fld, is_rot
+ real(kind = kreal), intent(in) :: r_ICB(0:2)
+ real(kind = kreal), intent(in) :: coef_fdm_free_ICB_vp2(0:1,3)
+ real(kind = kreal), intent(in) :: coef_fdm_free_ICB_vt2(0:1,3)
!
!
real(kind = kreal) :: d2s_dr2, d1t_dr1
@@ -108,10 +139,10 @@
!
!
!$omp parallel do private(inod,i_p1,i_p2,j,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
!
d2s_dr2 = coef_fdm_free_ICB_vp2( 0,3) * d_rj(inod,is_fld ) &
& + coef_fdm_free_ICB_vp2( 1,3) * d_rj(i_p1,is_fld )
@@ -120,8 +151,7 @@
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)*r_ICB(2) * d_rj(inod,is_fld )
end do
!$omp end parallel do
!
@@ -129,21 +159,25 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_free_diffuse2(coef_d, &
- & is_fld, is_diffuse)
+ subroutine cal_sph_nod_icb_free_diffuse2(jmax, kr_in, r_ICB, &
+ & coef_fdm_free_ICB_vp2, coef_fdm_free_ICB_vt2, &
+ & coef_d, is_fld, is_diffuse)
!
- 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) :: r_ICB(0:2)
+ real(kind = kreal), intent(in) :: coef_fdm_free_ICB_vp2(0:1,3)
+ real(kind = kreal), intent(in) :: coef_fdm_free_ICB_vt2(0:1,3)
!
real(kind = kreal) :: d2s_dr2, d2t_dr2
integer(kind = kint) :: inod, j, i_p1
!
!
!$omp parallel do private(inod,i_p1,j,d2s_dr2,d2t_dr2)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_ICB-1) * nidx_rj(2)
- i_p1 = inod + nidx_rj(2)
+ do j = 1, jmax
+ inod = j + (kr_in-1) * jmax
+ i_p1 = inod + jmax
!
d2s_dr2 = coef_fdm_free_ICB_vp2( 0,3) * d_rj(inod,is_fld ) &
& + coef_fdm_free_ICB_vp2( 1,3) * d_rj(i_p1,is_fld )
@@ -151,11 +185,9 @@
& + coef_fdm_free_ICB_vt2( 1,3) * d_rj(i_p1,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)*r_ICB(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)*r_ICB(2) * d_rj(inod,is_fld+2) )
end do
!$omp end parallel do
!
@@ -163,37 +195,37 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_free_w_diffuse2(coef_d, &
- & is_fld, is_diffuse)
+ subroutine cal_sph_nod_icb_free_w_diffuse2(jmax, kr_in, r_ICB, &
+ & fdm2_fix_fld_ICB, coef_fdm_free_ICB_vt2, &
+ & 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) :: r_ICB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_ICB(0:2,3)
+ real(kind = kreal), intent(in) :: coef_fdm_free_ICB_vt2(0:1,3)
!
integer(kind = kint) :: inod, j, i_p1, i_p2
real(kind = kreal) :: d2s_dr2,d2t_dr2
!
!
!$omp parallel do private(inod,i_p1,i_p2,j,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)
+ do j = 1, jmax
+ inod = j + (kr_in-1) * jmax
+ i_p1 = inod + jmax
+ i_p2 = i_p1 + jmax
!
d2s_dr2 = coef_fdm_free_ICB_vt2( 0,3) * d_rj(inod,is_fld ) &
& + coef_fdm_free_ICB_vt2( 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)
+ 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)*r_ICB(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)*r_ICB(2) * d_rj(inod,is_fld+2) )
end do
!$omp end parallel do
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_ICB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_ICB.f90
index ef1c356..733435a 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_ICB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_ICB.f90
@@ -20,6 +20,18 @@
!! & fdm2_fix_fld_ICB, coef_d, is_fld, is_diffuse)
!!@endverbatim
!!
+!!@n @param idx_rj_degree_zero Local address for degree 0
+!!@n @param idx_rj_degree_one(-1:1) Local address for degree 1
+!!@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 r_ICB(0:2) Radius at ICB
+!!@n @param Vt_ICB(jmax) Spectr data for toroidal velocity ICB
+!!
+!!@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 Address of poloidal velocity in d_rj
!!@n @param is_rot Address of poloidal vorticity in d_rj
More information about the CIG-COMMITS
mailing list