[cig-commits] [commit] Hiro_latest: Use boundary structures to apply perturbation of temperature (0d43a43)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Jan 28 11:13:31 PST 2014
Repository : ssh://geoshell/calypso
On branch : Hiro_latest
Link : https://github.com/geodynamics/calypso/compare/3c548304673360ddedd7d68c8095b3fb74a2b9ce...c4a1173efe5e196ae6c9f3494a66a74b5e60569a
>---------------------------------------------------------------
commit 0d43a43bae528ecc891c9cb890ef5e76ea93b03c
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Wed Nov 20 11:43:43 2013 -0800
Use boundary structures to apply perturbation of temperature
>---------------------------------------------------------------
0d43a43bae528ecc891c9cb890ef5e76ea93b03c
.../MHD_src/sph_MHD/cal_explicit_terms.f90 | 51 ++++++++-----
.../MHD_src/sph_MHD/cal_momentum_eq_explicit.f90 | 24 ++++---
.../MHD_src/sph_MHD/cal_nonlinear.f90 | 8 ++-
.../MHD_src/sph_MHD/cal_nonlinear_sph_MHD.f90 | 13 ++--
.../MHD_src/sph_MHD/cal_sol_sph_fluid_crank.f90 | 14 ++--
.../MHD_src/sph_MHD/cal_vorticity_terms_adams.f90 | 26 ++++---
.../MHD_src/sph_MHD/const_coriolis_sph.f90 | 2 +-
.../MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90 | 4 +-
.../MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90 | 19 +++--
.../MHD_src/sph_MHD/m_coef_fdm_free_CMB.f90 | 12 ++--
.../MHD_src/sph_MHD/m_coef_fdm_free_ICB.f90 | 12 ++--
.../MHD_src/sph_MHD/m_radial_matrices_sph.f90 | 10 ---
.../MHD_src/sph_MHD/set_bc_sph_mhd.f90 | 12 ++--
.../MHD_src/sph_MHD/set_radial_mat_sph.f90 | 36 ++++++++++
.../MHD_src/sph_MHD/set_reference_sph_mhd.f90 | 83 ++++++++++++----------
.../MHD_src/sph_MHD/trans_sph_velo_4_coriolis.f90 | 18 +++--
16 files changed, 220 insertions(+), 124 deletions(-)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_explicit_terms.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_explicit_terms.f90
index c1469df..db7edff 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_explicit_terms.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_explicit_terms.f90
@@ -11,16 +11,19 @@
!! subroutine cal_diff_induction_wSGS_adams
!! subroutine cal_diff_induction_MHD_euler
!!
-!! subroutine sel_heat_diff_adv_src_adams
-!! subroutine sel_heat_diff_adv_src_euler
+!! subroutine sel_heat_diff_adv_src_adams(kr_in, kr_out)
+!! subroutine sel_heat_diff_adv_src_euler(kr_in, kr_out)
!!
-!! subroutine sel_light_diff_adv_src_adams
-!! subroutine sel_light_diff_adv_src_euler
+!! subroutine sel_light_diff_adv_src_adams(kr_in, kr_out)
+!! subroutine sel_light_diff_adv_src_euler(kr_in, kr_out)
!!
!! subroutine set_ini_adams_mag_induct
-!! subroutine sel_ini_adams_heat_w_src
-!! subroutine sel_ini_adams_light_w_src
+!! subroutine sel_ini_adams_heat_w_src(kr_in, kr_out)
+!! subroutine sel_ini_adams_light_w_src(kr_in, kr_out)
!!@endverbatim
+!!
+!!@n @param kr_in Radial ID for inner boundary
+!!@n @param kr_out Radial ID for outer boundary
!
module cal_explicit_terms
!
@@ -137,13 +140,15 @@
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
!
- subroutine sel_heat_diff_adv_src_adams
+ subroutine sel_heat_diff_adv_src_adams(kr_in, kr_out)
!
use m_physical_property
use cal_diff_adv_src_explicit
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
!
- call sel_scalar_diff_adv_src_adams(nlayer_ICB, nlayer_CMB, &
+ call sel_scalar_diff_adv_src_adams(kr_in, kr_out, &
& ipol%i_t_diffuse, ipol%i_h_advect, ipol%i_heat_source, &
& ipol%i_temp, ipol%i_pre_heat, coef_exp_t, coef_h_src)
!
@@ -151,13 +156,15 @@
!
! ----------------------------------------------------------------------
!
- subroutine sel_heat_diff_adv_src_euler
+ subroutine sel_heat_diff_adv_src_euler(kr_in, kr_out)
!
use m_physical_property
use cal_diff_adv_src_explicit
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
!
- call sel_scalar_diff_adv_src_euler(nlayer_ICB, nlayer_CMB, &
+!
+ call sel_scalar_diff_adv_src_euler(kr_in, kr_out, &
& ipol%i_t_diffuse, ipol%i_h_advect, ipol%i_heat_source, &
& ipol%i_temp, coef_exp_t, coef_h_src)
!
@@ -166,13 +173,15 @@
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
!
- subroutine sel_light_diff_adv_src_adams
+ subroutine sel_light_diff_adv_src_adams(kr_in, kr_out)
!
use m_physical_property
use cal_diff_adv_src_explicit
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
!
- call sel_scalar_diff_adv_src_adams(nlayer_ICB, nlayer_CMB, &
+ call sel_scalar_diff_adv_src_adams(kr_in, kr_out, &
& ipol%i_c_diffuse, ipol%i_c_advect, ipol%i_light_source, &
& ipol%i_light, ipol%i_pre_composit, coef_exp_c, coef_c_src)
!
@@ -180,13 +189,15 @@
!
! ----------------------------------------------------------------------
!
- subroutine sel_light_diff_adv_src_euler
+ subroutine sel_light_diff_adv_src_euler(kr_in, kr_out)
!
use m_physical_property
use cal_diff_adv_src_explicit
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
!
- call sel_scalar_diff_adv_src_euler(nlayer_ICB, nlayer_CMB, &
+!
+ call sel_scalar_diff_adv_src_euler(kr_in, kr_out, &
& ipol%i_c_diffuse, ipol%i_c_advect, ipol%i_light_source, &
& ipol%i_light, coef_exp_c, coef_c_src)
!
@@ -211,13 +222,15 @@
!
! ----------------------------------------------------------------------
!
- subroutine sel_ini_adams_heat_w_src
+ subroutine sel_ini_adams_heat_w_src(kr_in, kr_out)
!
use m_physical_property
use cal_diff_adv_src_explicit
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
!
- call sel_ini_adams_scalar_w_src(nlayer_ICB, nlayer_CMB, &
+ call sel_ini_adams_scalar_w_src(kr_in, kr_out, &
& ipol%i_h_advect, ipol%i_heat_source, ipol%i_pre_heat, &
& coef_h_src)
!
@@ -225,13 +238,15 @@
!
! ----------------------------------------------------------------------
!
- subroutine sel_ini_adams_light_w_src
+ subroutine sel_ini_adams_light_w_src(kr_in, kr_out)
!
use m_physical_property
use cal_diff_adv_src_explicit
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
!
- call sel_ini_adams_scalar_w_src(nlayer_ICB, nlayer_CMB, &
+ call sel_ini_adams_scalar_w_src(kr_in, kr_out, &
& ipol%i_c_advect, ipol%i_light_source, ipol%i_pre_composit, &
& coef_c_src)
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_momentum_eq_explicit.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_momentum_eq_explicit.f90
index 3ad85b8..1832d24 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_momentum_eq_explicit.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_momentum_eq_explicit.f90
@@ -53,6 +53,7 @@
!
subroutine cal_expricit_sph_adams
!
+ use m_boundary_params_sph_MHD
use cal_explicit_terms
use cal_vorticity_terms_adams
use cal_nonlinear_sph_MHD
@@ -60,17 +61,19 @@
!
!$omp parallel
if(iflag_t_evo_4_velo .gt. id_no_evolution) then
- call cal_vorticity_eq_adams
+ call cal_vorticity_eq_adams(sph_bc_U%kr_in, sph_bc_U%kr_out)
end if
!
if(iflag_t_evo_4_magne .gt. id_no_evolution) then
call cal_diff_induction_MHD_adams
end if
if(iflag_t_evo_4_temp .gt. id_no_evolution) then
- call sel_heat_diff_adv_src_adams
+ call sel_heat_diff_adv_src_adams &
+ & (sph_bc_T%kr_in, sph_bc_T%kr_out)
end if
if(iflag_t_evo_4_composit .gt. id_no_evolution) then
- call sel_light_diff_adv_src_adams
+ call sel_light_diff_adv_src_adams &
+ & (sph_bc_C%kr_in, sph_bc_C%kr_out)
end if
!$omp end parallel
!
@@ -80,6 +83,7 @@
!
subroutine cal_expricit_sph_euler(i_step)
!
+ use m_boundary_params_sph_MHD
use cal_explicit_terms
use cal_vorticity_terms_adams
!
@@ -88,17 +92,19 @@
!
!$omp parallel
if(iflag_t_evo_4_velo .gt. id_no_evolution) then
- call cal_vorticity_eq_euler
+ call cal_vorticity_eq_euler(sph_bc_U%kr_in, sph_bc_U%kr_out)
end if
!
if(iflag_t_evo_4_temp .gt. id_no_evolution) then
- call sel_heat_diff_adv_src_euler
+ call sel_heat_diff_adv_src_euler &
+ & (sph_bc_T%kr_in, sph_bc_T%kr_out)
end if
if(iflag_t_evo_4_magne .gt. id_no_evolution) then
call cal_diff_induction_MHD_euler
end if
if(iflag_t_evo_4_composit .gt. id_no_evolution) then
- call sel_light_diff_adv_src_euler
+ call sel_light_diff_adv_src_euler &
+ & (sph_bc_C%kr_in, sph_bc_C%kr_out)
end if
!
if (i_step .eq. 1) then
@@ -106,13 +112,15 @@
call set_ini_adams_inertia
end if
if(iflag_t_evo_4_temp .gt. id_no_evolution) then
- call sel_ini_adams_heat_w_src
+ call sel_ini_adams_heat_w_src &
+ & (sph_bc_T%kr_in, sph_bc_T%kr_out)
end if
if(iflag_t_evo_4_magne .gt. id_no_evolution) then
call set_ini_adams_mag_induct
end if
if(iflag_t_evo_4_composit .gt. id_no_evolution) then
- call sel_ini_adams_light_w_src
+ call sel_ini_adams_light_w_src &
+ & (sph_bc_C%kr_in, sph_bc_C%kr_out)
end if
end if
!$omp end parallel
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_nonlinear.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_nonlinear.f90
index 32b2b72..d11f565 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_nonlinear.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_nonlinear.f90
@@ -33,6 +33,7 @@
subroutine nonlinear
!
use m_sph_phys_address
+ use m_boundary_params_sph_MHD
use cal_vorticity_terms_adams
use const_coriolis_sph
!
@@ -47,7 +48,8 @@
call nonlinear_by_pseudo_sph
!
if (iflag_4_ref_temp .eq. id_sphere_ref_temp) then
- call add_reftemp_advect_sph_MHD
+ call add_reftemp_advect_sph_MHD &
+ & (sph_bc_T%kr_in, sph_bc_T%kr_out)
end if
!
!* ---- set coriolis term
@@ -137,6 +139,7 @@
subroutine licv_exp
!
use m_sph_phys_address
+ use m_boundary_params_sph_MHD
use cal_nonlinear_sph_MHD
use cal_vorticity_terms_adams
use const_coriolis_sph
@@ -159,7 +162,8 @@
!$omp end parallel do
!
if (iflag_4_ref_temp .eq. id_sphere_ref_temp) then
- call add_reftemp_advect_sph_MHD
+ call add_reftemp_advect_sph_MHD &
+ & (sph_bc_T%kr_in, sph_bc_T%kr_out)
end if
!
!$omp parallel
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_nonlinear_sph_MHD.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_nonlinear_sph_MHD.f90
index 87aa60f..56c35e5 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_nonlinear_sph_MHD.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_nonlinear_sph_MHD.f90
@@ -8,8 +8,11 @@
!!
!!@verbatim
!! subroutine s_cal_nonlinear_sph_MHD
-!! subroutine add_reftemp_advect_sph_MHD
+!! subroutine add_reftemp_advect_sph_MHD(kr_in, kr_out)
!!@endverbatim
+!!
+!!@n @param kr_in Radial ID for inner boundary
+!!@n @param kr_out Radial ID for outer boundary
!
module cal_nonlinear_sph_MHD
!
@@ -77,15 +80,17 @@
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
- subroutine add_reftemp_advect_sph_MHD
+ subroutine add_reftemp_advect_sph_MHD(kr_in, kr_out)
!
use m_schmidt_poly_on_rtm
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
integer(kind= kint) :: ist, ied, inod, j, k
!
!
- ist = (nlayer_ICB-1)*nidx_rj(2) + 1
- ied = nlayer_CMB * nidx_rj(2)
+ ist = (kr_in-1)*nidx_rj(2) + 1
+ ied = kr_out * nidx_rj(2)
!$omp parallel do private (inod,j,k)
do inod = ist, ied
j = mod((inod-1),nidx_rj(2)) + 1
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 52a50dc..90e929c 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
@@ -159,7 +159,7 @@
& nidx_rj(2), nidx_rj(1), p_poisson_lu, i_p_pivot, &
& d_rj(1,ipol%i_press) )
!
- call adjust_by_ave_pressure_on_CMB
+ call adjust_by_ave_pressure_on_CMB(sph_bc_U%kr_out)
!
! if (idx_rj_degree_zero .gt. 0) then
! icmb = (sph_bc_U%kr_out-1)*nidx_rj(2) + idx_rj_degree_zero
@@ -247,7 +247,7 @@
& sph_bc_T%ICB_flux, coef_d_temp, coef_imp_t, dt, &
& ipol%i_temp)
else
- call set_fixed_scalar_sph(nidx_rj(2), ione, nlayer_ICB, &
+ call set_fixed_scalar_sph(nidx_rj(2), ione, sph_bc_T%kr_in, &
& ipol%i_temp, sph_bc_T%ICB_fld)
end if
!
@@ -257,8 +257,9 @@
& sph_bc_T%CMB_flux, coef_d_temp, coef_imp_t, dt, &
& ipol%i_temp)
else
- call set_fixed_scalar_sph(nidx_rj(2), nlayer_CMB, nidx_rj(1), &
- & ipol%i_temp, sph_bc_T%CMB_fld)
+ call set_fixed_scalar_sph &
+ & (nidx_rj(2), sph_bc_T%kr_out, nidx_rj(1), &
+ & ipol%i_temp, sph_bc_T%CMB_fld)
end if
!
call lubksb_3band_mul(np_smp, idx_rj_smp_stack(0,2), &
@@ -283,7 +284,7 @@
& sph_bc_C%ICB_flux, coef_d_light, coef_imp_c, dt, &
& ipol%i_light)
else
- call set_fixed_scalar_sph(nidx_rj(2), ione, nlayer_ICB, &
+ call set_fixed_scalar_sph(nidx_rj(2), ione, sph_bc_C%kr_in, &
& ipol%i_light, sph_bc_C%ICB_fld)
end if
!
@@ -293,7 +294,8 @@
& sph_bc_C%CMB_flux, coef_d_light, coef_imp_c, dt, &
& ipol%i_light)
else
- call set_fixed_scalar_sph(nidx_rj(2), nlayer_CMB, nidx_rj(1), &
+ call set_fixed_scalar_sph &
+ & (nidx_rj(2), sph_bc_C%kr_out, nidx_rj(1), &
& ipol%i_light, sph_bc_C%CMB_fld)
end if
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_vorticity_terms_adams.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_vorticity_terms_adams.f90
index 66010cb..ee5a472 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_vorticity_terms_adams.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_vorticity_terms_adams.f90
@@ -7,8 +7,8 @@
!>@brief Evoluve the vorticity equation by explicit scheme
!!
!!@verbatim
-!! subroutine cal_vorticity_eq_adams
-!! subroutine cal_vorticity_eq_euler
+!! subroutine cal_vorticity_eq_adams(kr_in, kr_out)
+!! subroutine cal_vorticity_eq_euler(kr_in, kr_out)
!!
!! subroutine set_MHD_terms_to_force(it_rot_buo)
!! subroutine set_rot_cv_terms_to_force(it_rot_buo)
@@ -20,6 +20,12 @@
!!
!! subroutine set_ini_adams_inertia
!!@endverbatim
+!!
+!!@n @param kr_in Radial ID for inner boundary
+!!@n @param kr_out Radial ID for outer boundary
+!!
+!!@n @param it_rot_buo Spectr field address
+!! for toroidal curl of buodyancy
!
module cal_vorticity_terms_adams
!
@@ -38,15 +44,17 @@
!
! ----------------------------------------------------------------------
!
- subroutine cal_vorticity_eq_adams
+ subroutine cal_vorticity_eq_adams(kr_in, kr_out)
!
use m_physical_property
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
integer(kind = kint) :: inod, ist, ied
!
!
- ist = (nlayer_ICB-1)*nidx_rj(2) + 1
- ied = nlayer_CMB * nidx_rj(2)
+ ist = (kr_in-1)*nidx_rj(2) + 1
+ ied = kr_out * nidx_rj(2)
!$omp do private (inod)
do inod = ist, ied
d_rj(inod,ipol%i_vort) = d_rj(inod,ipol%i_vort) &
@@ -67,15 +75,17 @@
!
! ----------------------------------------------------------------------
!
- subroutine cal_vorticity_eq_euler
+ subroutine cal_vorticity_eq_euler(kr_in, kr_out)
!
use m_physical_property
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
integer(kind = kint) :: inod, ist, ied
!
!
- ist = (nlayer_ICB-1)*nidx_rj(2) + 1
- ied = nlayer_CMB * nidx_rj(2)
+ ist = (kr_in-1)*nidx_rj(2) + 1
+ ied = kr_out * nidx_rj(2)
!$omp do private (inod)
do inod = ist, ied
d_rj(inod,ipol%i_vort) = d_rj(inod,ipol%i_vort) &
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_coriolis_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_coriolis_sph.f90
index b4ad5f8..7d8b236 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_coriolis_sph.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_coriolis_sph.f90
@@ -81,7 +81,7 @@
!
!
if (iflag_debug.eq.1) write(*,*) 's_trans_sph_velo_4_coriolis'
- call s_trans_sph_velo_4_coriolis
+ call s_trans_sph_velo_4_coriolis(sph_bc_U%kr_in, sph_bc_U%kr_out)
!
call s_sum_rot_coriolis_rj_sph(sph_bc_U%kr_in, sph_bc_U%kr_out, &
& coef_cor)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90
index b1031e1..089b8a4 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90
@@ -43,7 +43,7 @@
!
!$omp parallel
call set_radial_scalar_evo_mat_sph(nidx_rj(1), nidx_rj(2), &
- & nlayer_ICB, nlayer_CMB, coef_imp_t, coef_d_temp, &
+ & sph_bc_T%kr_in, sph_bc_T%kr_out, coef_imp_t, coef_d_temp, &
& temp_evo_mat)
!$omp end parallel
!
@@ -91,7 +91,7 @@
!
!$omp parallel
call set_radial_scalar_evo_mat_sph(nidx_rj(1), nidx_rj(2), &
- & nlayer_ICB, nlayer_CMB, coef_imp_c, coef_d_light, &
+ & sph_bc_C%kr_in, sph_bc_C%kr_out, coef_imp_c, coef_d_light, &
& composit_evo_mat)
!$omp end parallel
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90
index 0e00b82..553a7a1 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90
@@ -48,14 +48,23 @@
!
!
!$omp parallel
+ call set_unit_umat_4_poisson(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%kr_out, vs_poisson_mat)
+ call set_unit_umat_4_poisson(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%kr_out, p_poisson_mat)
+!$omp end parallel
+!
+!$omp parallel
call set_radial_vect_evo_mat_sph(nidx_rj(1), nidx_rj(2), &
- & nlayer_ICB, nlayer_CMB, coef_imp_v, coef_d_velo, vt_evo_mat)
+ & sph_bc_U%kr_in, sph_bc_U%kr_out, coef_imp_v, coef_d_velo, &
+ & vt_evo_mat)
call set_radial_vect_evo_mat_sph(nidx_rj(1), nidx_rj(2), &
- & nlayer_ICB, nlayer_CMB, coef_imp_v, coef_d_velo, wt_evo_mat)
+ & sph_bc_U%kr_in, sph_bc_U%kr_out, coef_imp_v, coef_d_velo, &
+ & wt_evo_mat)
call set_radial_vp3_mat_sph(nidx_rj(1), nidx_rj(2), &
- & nlayer_ICB, nlayer_CMB, vs_poisson_mat)
+ & sph_bc_U%kr_in, sph_bc_U%kr_out, vs_poisson_mat)
call set_radial_press_mat_sph(nidx_rj(1), nidx_rj(2), &
- & nlayer_ICB, nlayer_CMB, coef_press, p_poisson_mat)
+ & sph_bc_U%kr_in, sph_bc_U%kr_out, coef_press, p_poisson_mat)
!$omp end parallel
!
! Boundary condition for ICB
@@ -114,7 +123,7 @@
!
!
call cal_mat_product_3band_mul(nidx_rj(1), nidx_rj(2), &
- & nlayer_ICB, nlayer_CMB, wt_evo_mat, vs_poisson_mat, &
+ & sph_bc_U%kr_in, sph_bc_U%kr_out, wt_evo_mat, vs_poisson_mat, &
& vp_evo_mat)
!
if(i_debug .eq. iflag_full_msg) &
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_free_CMB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_free_CMB.f90
index 4c75544..83dfcdf 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_free_CMB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_free_CMB.f90
@@ -8,8 +8,8 @@
!! at CMB with free slip boundary
!!
!!@verbatim
-!! subroutine cal_2nd_CMB_free_vp_bc_fdm(r_from_CMB)
-!! subroutine cal_2nd_CMB_free_vt_bc_fdm(r_from_CMB)
+!! subroutine cal_fdm2_CMB_free_vp(r_from_CMB)
+!! subroutine cal_fdm2_CMB_free_vt(r_from_CMB)
!!
!! subroutine check_coef_fdm_free_CMB
!!
@@ -82,7 +82,7 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_2nd_CMB_free_vp_bc_fdm(r_from_CMB)
+ subroutine cal_fdm2_CMB_free_vp(r_from_CMB)
!
real(kind = kreal) :: r_from_CMB(-1:0)
!
@@ -122,11 +122,11 @@
fdm2_free_vp_CMB(0, 3) = mat_fdm_CMB_free_vp(3,1)
fdm2_free_vp_CMB(-1,3) = mat_fdm_CMB_free_vp(3,3)
!
- end subroutine cal_2nd_CMB_free_vp_bc_fdm
+ end subroutine cal_fdm2_CMB_free_vp
!
! -----------------------------------------------------------------------
!
- subroutine cal_2nd_CMB_free_vt_bc_fdm(r_from_CMB)
+ subroutine cal_fdm2_CMB_free_vt(r_from_CMB)
!
real(kind = kreal) :: r_from_CMB(-1:0)
!
@@ -166,7 +166,7 @@
fdm2_free_vt_CMB(0, 3) = mat_fdm_CMB_free_vt(3,1)
fdm2_free_vt_CMB(-1,3) = mat_fdm_CMB_free_vt(3,3)
!
- end subroutine cal_2nd_CMB_free_vt_bc_fdm
+ end subroutine cal_fdm2_CMB_free_vt
!
! -----------------------------------------------------------------------
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_free_ICB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_free_ICB.f90
index 06ed07a..d1b320d 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_free_ICB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_free_ICB.f90
@@ -8,8 +8,8 @@
!! at CMB with free slip boundary
!!
!!@verbatim
-!! subroutine cal_2nd_ICB_free_vp_bc_fdm(r_from_ICB)
-!! subroutine cal_2nd_ICB_free_vt_bc_fdm(r_from_ICB)
+!! subroutine cal_fdm2_ICB_free_vp(r_from_ICB)
+!! subroutine cal_fdm2_ICB_free_vt(r_from_ICB)
!!
!! subroutine check_coef_fdm_free_ICB
!!
@@ -82,7 +82,7 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_2nd_ICB_free_vp_bc_fdm(r_from_ICB)
+ subroutine cal_fdm2_ICB_free_vp(r_from_ICB)
!
real(kind = kreal), intent(in) :: r_from_ICB(0:1)
!
@@ -121,11 +121,11 @@
fdm2_free_vp_ICB(0,3) = mat_fdm_ICB_free_vp(3,1)
fdm2_free_vp_ICB(1,3) = mat_fdm_ICB_free_vp(3,3)
!
- end subroutine cal_2nd_ICB_free_vp_bc_fdm
+ end subroutine cal_fdm2_ICB_free_vp
!
! -----------------------------------------------------------------------
!
- subroutine cal_2nd_ICB_free_vt_bc_fdm(r_from_ICB)
+ subroutine cal_fdm2_ICB_free_vt(r_from_ICB)
!
real(kind = kreal), intent(in) :: r_from_ICB(0:1)
!
@@ -164,7 +164,7 @@
fdm2_free_vt_ICB(0,3) = mat_fdm_ICB_free_vt(3,1)
fdm2_free_vt_ICB(1,3) = mat_fdm_ICB_free_vt(3,3)
!
- end subroutine cal_2nd_ICB_free_vt_bc_fdm
+ end subroutine cal_fdm2_ICB_free_vt
!
! -----------------------------------------------------------------------
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/m_radial_matrices_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/m_radial_matrices_sph.f90
index 873e2c1..9bb5645 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/m_radial_matrices_sph.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/m_radial_matrices_sph.f90
@@ -172,16 +172,6 @@
vt_evo_mat(2,1:nri,1:jmax) = 1.0d0
wt_evo_mat(2,1:nri,1:jmax) = 1.0d0
!
- if(nlayer_ICB .gt. 1) then
- p_poisson_mat(2,1:nlayer_ICB-1,1:jmax) = 1.0d0
- vs_poisson_mat(2,1:nlayer_ICB-1,1:jmax) = 1.0d0
- end if
-!
- if(nlayer_CMB .lt. nri) then
- p_poisson_mat(2,nlayer_CMB+1:nri,1:jmax) = 1.0d0
- vs_poisson_mat(2,nlayer_CMB+1:nri,1:jmax) = 1.0d0
- end if
-!
end subroutine allocate_velo_mat_sph
!
! -----------------------------------------------------------------------
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 f7f2385..7d24e61 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
@@ -75,11 +75,11 @@
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_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_fdm2_ICB_free_vp(radius_1d_rj_r(sph_bc_U%kr_in))
+ call cal_fdm2_ICB_free_vt(radius_1d_rj_r(sph_bc_U%kr_in))
!
- 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_fdm2_CMB_free_vp(radius_1d_rj_r(sph_bc_U%kr_out-1))
+ call cal_fdm2_CMB_free_vt(radius_1d_rj_r(sph_bc_U%kr_out-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))
@@ -88,8 +88,8 @@
! Set reference temperature and adjust boundary conditions
!
call allocate_reft_rj_data
- call set_ref_temp_sph_mhd
- call adjust_sph_temp_bc_by_reftemp
+ call set_ref_temp_sph_mhd(sph_bc_T)
+ call adjust_sph_temp_bc_by_reftemp(sph_bc_T)
!
! Check data
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_radial_mat_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_radial_mat_sph.f90
index d8bb450..9a548d4 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_radial_mat_sph.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_radial_mat_sph.f90
@@ -17,6 +17,9 @@
!! subroutine set_radial_press_mat_sph(nri, jmax, kr_in, kr_out, &
!! & coef_p, poisson_mat)
!!
+!! subroutine set_unit_umat_4_poisson(nri, jmax, &
+!! & kr_in, kr_out, poisson_mat)
+!!
!! Format of band matrix
!! | a(2,1) a(1,2) ........ 0 0 |
!! | a(3,1) a(2,2) ........ . . |
@@ -192,4 +195,37 @@
!
! -----------------------------------------------------------------------
!
+ subroutine set_unit_umat_4_poisson(nri, jmax, &
+ & kr_in, kr_out, poisson_mat)
+!
+ integer(kind = kint), intent(in) :: jmax, nri
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
+ real(kind = kreal), intent(inout) :: poisson_mat(3,nri,jmax)
+!
+ integer(kind = kint) :: k, j
+!
+!
+!$omp do private (k,j)
+ do j = 1, jmax
+ do k = 1, nri
+ poisson_mat(3,k,j) = zero
+ poisson_mat(1,k,j) = zero
+ end do
+ do k = 1, kr_in-1
+ poisson_mat(2,k,j) = one
+ end do
+ do k = kr_in, kr_out
+ poisson_mat(2,k,j) = zero
+ end do
+ do k = kr_out+1, nri
+ poisson_mat(2,k,j) = one
+ end do
+ end do
+!$omp end do nowait
+!
+ end subroutine set_unit_umat_4_poisson
+!
+! -----------------------------------------------------------------------
+!
end module set_radial_mat_sph
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_reference_sph_mhd.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_reference_sph_mhd.f90
index e46b913..8b6b759 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_reference_sph_mhd.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_reference_sph_mhd.f90
@@ -7,9 +7,10 @@
!> @brief Convert temperature data using reference temperature
!!
!!@verbatim
-!! subroutine adjust_by_ave_pressure_on_CMB
+!! subroutine adjust_by_ave_pressure_on_CMB(kr_out)
!!
-!! subroutine set_ref_temp_sph_mhd
+!! subroutine set_ref_temp_sph_mhd(sph_bc_T)
+!! subroutine adjust_sph_temp_bc_by_reftemp(sph_bc_T)
!!
!! subroutine sync_temp_by_per_temp_sph
!! d_rj(inod,ipol%i_temp): T => \Theta = T - T0
@@ -24,10 +25,11 @@
!! d_rj(inod,ipol%i_grad_t): d \Theta / dr => dT / dr
!! d_rj(inod,ipol%i_grad_part_t): d \Theta / dr
!!
-!! subroutine adjust_sph_temp_bc_by_reftemp
!! subroutine delete_zero_degree_comp(is_fld)
!!@endverbatim
!!
+!!@param sph_bc_T Structure for basic boundary condition parameters
+!! for temperature
!!@n @param is_fld Address of poloidal component
!
module set_reference_sph_mhd
@@ -48,17 +50,19 @@
!
! -----------------------------------------------------------------------
!
- subroutine adjust_by_ave_pressure_on_CMB
+ subroutine adjust_by_ave_pressure_on_CMB(kr_out)
!
use m_sph_phys_address
!
+ integer(kind = kint), intent(in) :: kr_out
+!
integer(kind = kint) :: k, inod
real(kind = kreal) :: ref_p
!
!
if (idx_rj_degree_zero .eq. 0) return
!
- inod = idx_rj_degree_zero + (nlayer_CMB-1)*nidx_rj(2)
+ inod = idx_rj_degree_zero + (kr_out-1)*nidx_rj(2)
ref_p = d_rj(inod,ipol%i_press)
!
do k = 1, nidx_rj(1)
@@ -71,18 +75,22 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine set_ref_temp_sph_mhd
+ subroutine set_ref_temp_sph_mhd(sph_bc_T)
+!
+ use t_boundary_params_sph_MHD
+!
+ type(sph_boundary_type), intent(in) :: sph_bc_T
!
integer (kind = kint) :: k
!
! set reference temperature (for spherical shell)
!
if (iflag_4_ref_temp .eq. id_sphere_ref_temp) then
- do k = 1, nlayer_ICB-1
+ do k = 1, sph_bc_T%kr_in-1
reftemp_rj(k,0) = high_temp
reftemp_rj(k,1) = zero
end do
- do k = nlayer_ICB, nlayer_CMB
+ do k = sph_bc_T%kr_in, sph_bc_T%kr_out
reftemp_rj(k,0) = (depth_high_t*depth_low_t*ar_1d_rj(k,1) &
& * (high_temp - low_temp) &
& - depth_high_t*high_temp &
@@ -92,7 +100,7 @@
& * (high_temp - low_temp) &
& / (depth_low_t - depth_high_t)
end do
- do k = nlayer_CMB+1, nidx_rj(1)
+ do k = sph_bc_T%kr_out+1, nidx_rj(1)
reftemp_rj(k,0) = low_temp
reftemp_rj(k,1) = zero
end do
@@ -107,6 +115,36 @@
!
! -----------------------------------------------------------------------
!
+ subroutine adjust_sph_temp_bc_by_reftemp(sph_bc_T)
+!
+ use m_spheric_parameter
+ use m_sph_spectr_data
+ use t_boundary_params_sph_MHD
+!
+ type(sph_boundary_type), intent(in) :: sph_bc_T
+!
+!
+ if(idx_rj_degree_zero .gt. 0 &
+ & .and. iflag_4_ref_temp .eq. id_sphere_ref_temp) then
+ sph_bc_T%ICB_fld(idx_rj_degree_zero) &
+ & = sph_bc_T%ICB_fld(idx_rj_degree_zero) &
+ & - reftemp_rj(sph_bc_T%kr_in,0)
+ sph_bc_T%CMB_fld(idx_rj_degree_zero) &
+ & = sph_bc_T%CMB_fld(idx_rj_degree_zero) &
+ & - reftemp_rj(sph_bc_T%kr_out,0)
+ sph_bc_T%ICB_flux(idx_rj_degree_zero) &
+ & = sph_bc_T%ICB_flux(idx_rj_degree_zero) &
+ & - reftemp_rj(sph_bc_T%kr_in,1)
+ sph_bc_T%CMB_flux(idx_rj_degree_zero) &
+ & = sph_bc_T%CMB_flux(idx_rj_degree_zero) &
+ & - reftemp_rj(sph_bc_T%kr_out,1)
+ end if
+!
+ end subroutine adjust_sph_temp_bc_by_reftemp
+!
+! -----------------------------------------------------------------------
+! -----------------------------------------------------------------------
+!
subroutine sync_temp_by_per_temp_sph
!
use m_sph_phys_address
@@ -168,33 +206,6 @@
end subroutine trans_per_temp_to_temp_sph
!
! -----------------------------------------------------------------------
-!
- subroutine adjust_sph_temp_bc_by_reftemp
-!
- use m_spheric_parameter
- use m_sph_spectr_data
- use m_boundary_params_sph_MHD
-!
-!
- if(idx_rj_degree_zero .gt. 0 &
- & .and. iflag_4_ref_temp .eq. id_sphere_ref_temp) then
- sph_bc_T%ICB_fld(idx_rj_degree_zero) &
- & = sph_bc_T%ICB_fld(idx_rj_degree_zero) &
- & - reftemp_rj(nlayer_ICB,0)
- sph_bc_T%CMB_fld(idx_rj_degree_zero) &
- & = sph_bc_T%CMB_fld(idx_rj_degree_zero) &
- & - reftemp_rj(nlayer_CMB,0)
- sph_bc_T%ICB_flux(idx_rj_degree_zero) &
- & = sph_bc_T%ICB_flux(idx_rj_degree_zero) &
- & - reftemp_rj(nlayer_ICB,1)
- sph_bc_T%CMB_flux(idx_rj_degree_zero) &
- & = sph_bc_T%CMB_flux(idx_rj_degree_zero) &
- & - reftemp_rj(nlayer_CMB,1)
- end if
-!
- end subroutine adjust_sph_temp_bc_by_reftemp
-!
-! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
subroutine delete_zero_degree_comp(is_fld)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/trans_sph_velo_4_coriolis.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/trans_sph_velo_4_coriolis.f90
index bee568e..2795ae1 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/trans_sph_velo_4_coriolis.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/trans_sph_velo_4_coriolis.f90
@@ -7,8 +7,11 @@
!>@brief Data transfer to evaluate Coriolis term
!!
!!@verbatim
-!! subroutine s_trans_sph_velo_4_coriolis
+!! subroutine s_trans_sph_velo_4_coriolis(kr_in, kr_out)
!!@endverbatim
+!!
+!!@n @param kr_in Radial ID for inner boundary
+!!@n @param kr_out Radial ID for outer boundary
!
module trans_sph_velo_4_coriolis
!
@@ -30,12 +33,14 @@
!
! -----------------------------------------------------------------------
!
- subroutine s_trans_sph_velo_4_coriolis
+ subroutine s_trans_sph_velo_4_coriolis(kr_in, kr_out)
!
use solver_sph_coriolis_sr
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
+!
!
- call copy_sph_diff_velocties
+ call copy_sph_diff_velocties(kr_in, kr_out)
!
call solver_sph_coriolis_sr_5(nshift_j_cor, &
& nidx_rj(2), nidx_rj(1), d_rj_cor(1,1), &
@@ -46,15 +51,16 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine copy_sph_diff_velocties
+ subroutine copy_sph_diff_velocties(kr_in, kr_out)
!
use m_schmidt_poly_on_rtm
!
+ integer(kind = kint), intent(in) :: kr_in, kr_out
integer(kind = kint) :: ist, ied, i, k, j
!
!
- ist = (nlayer_ICB-1)*nidx_rj(2) + 1
- ied = nlayer_CMB * nidx_rj(2)
+ ist = (kr_in-1)*nidx_rj(2) + 1
+ ied = kr_out * nidx_rj(2)
!$omp parallel do private(k,j)
do i = ist, ied
j = mod((i-1),nidx_rj(2)) + 1
More information about the CIG-COMMITS
mailing list