[cig-commits] [commit] Hiro_latest: Use structure for non-slip boundary for CMB (ee5fcce)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Nov 18 16:22:01 PST 2013
Repository : ssh://geoshell/calypso
On branch : Hiro_latest
Link : https://github.com/geodynamics/calypso/compare/93e9f8f974c7a247c8f02e54ec18de063f86c8fb...3c548304673360ddedd7d68c8095b3fb74a2b9ce
>---------------------------------------------------------------
commit ee5fcce818c290760fe056ca2edd3d61a566826e
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Thu Nov 14 19:17:28 2013 -0800
Use structure for non-slip boundary for CMB
>---------------------------------------------------------------
ee5fcce818c290760fe056ca2edd3d61a566826e
.../MHD_src/sph_MHD/cal_sol_sph_fluid_crank.f90 | 10 +-
.../MHD_src/sph_MHD/const_sph_diffusion.f90 | 22 ++-
.../MHD_src/sph_MHD/const_sph_radial_grad.f90 | 32 +++-
.../MHD_src/sph_MHD/const_sph_rotation.f90 | 37 +++--
.../MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90 | 162 ++++++++++---------
.../MHD_src/sph_MHD/set_sph_exp_rigid_ICB.f90 | 174 +++++++++++----------
6 files changed, 259 insertions(+), 178 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 e478b73..e5af4f7 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
@@ -75,15 +75,19 @@
if (sph_bc_U%iflag_icb .eq. iflag_free_slip) then
call cal_sph_nod_icb_free_vpol2(ipol%i_velo)
else if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then
- call cal_sph_nod_icb_rotate_velo2(ipol%i_velo)
+ call cal_sph_nod_icb_rotate_velo2 &
+ & (idx_rj_degree_zero, idx_rj_degree_one, nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%r_ICB, vt_ICB_bc, ipol%i_velo)
else
- call cal_sph_nod_icb_rigid_velo2(ipol%i_velo)
+ call cal_sph_nod_icb_rigid_velo2(nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%r_ICB, vt_ICB_bc, ipol%i_velo)
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
call cal_sph_nod_cmb_free_vpol2(ipol%i_velo)
else
- call cal_sph_nod_cmb_rigid_velo2(ipol%i_velo)
+ 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)
end if
!
! write(my_rank+70,*) 'k, j, inod, vp_rhs, vt_rhs'
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 ce758db..d88c010 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
@@ -69,8 +69,10 @@
call cal_sph_nod_icb_free_diffuse2(coef_d_velo, &
& ipol%i_velo, ipol%i_v_diffuse)
else
- call cal_sph_nod_icb_rigid_diffuse2(coef_d_velo, &
- & ipol%i_velo, ipol%i_v_diffuse)
+ call cal_sph_nod_icb_rigid_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & sph_bc_U%fdm2_fix_fld_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & coef_d_velo, ipol%i_velo, ipol%i_v_diffuse)
end if
call cal_dsdr_sph_no_bc_in_2(nidx_rj(2), &
& sph_bc_U%kr_in, sph_bc_U%fdm2_fix_fld_ICB, &
@@ -86,8 +88,10 @@
call cal_sph_nod_cmb_free_diffuse2(coef_d_velo, &
& ipol%i_velo, ipol%i_v_diffuse)
else
- call cal_sph_nod_cmb_rigid_diffuse2(coef_d_velo, &
- & ipol%i_velo, ipol%i_v_diffuse)
+ call cal_sph_nod_cmb_rigid_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & sph_bc_U%fdm2_fix_fld_CMB, sph_bc_U%fdm2_fix_dr_CMB, &
+ & coef_d_velo, ipol%i_velo, ipol%i_v_diffuse)
end if
call cal_dsdr_sph_no_bc_out_2(nidx_rj(2), &
& sph_bc_U%kr_out, sph_bc_U%fdm2_fix_fld_CMB, &
@@ -123,8 +127,9 @@
call cal_sph_nod_icb_free_w_diffuse2(coef_d_velo, &
& ipol%i_vort, ipol%i_w_diffuse)
else
- call cal_sph_nod_icb_rgd_w_diffuse2(coef_d_velo, &
- & ipol%i_vort, ipol%i_w_diffuse)
+ 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, &
+ & coef_d_velo, ipol%i_vort, ipol%i_w_diffuse)
end if
!
if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then
@@ -140,8 +145,9 @@
call cal_sph_nod_cmb_free_w_diffuse2(coef_d_velo, &
& ipol%i_vort, ipol%i_w_diffuse)
else
- call cal_sph_nod_cmb_rgd_w_diffuse2(coef_d_velo, &
- & ipol%i_vort, ipol%i_w_diffuse)
+ 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, &
+ & coef_d_velo, ipol%i_vort, ipol%i_w_diffuse)
end if
!
call cal_dsdr_sph_no_bc_out_2(nidx_rj(2), &
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 3208bc2..605a871 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
@@ -140,17 +140,29 @@
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)
else if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then
- call cal_sph_nod_icb_rotate_velo2(ipol%i_velo)
- call cal_sph_nod_icb_rigid_rot2(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_icb_rotate_velo2 &
+ & (idx_rj_degree_zero, idx_rj_degree_one, nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%r_ICB, vt_ICB_bc, ipol%i_velo)
+ call cal_sph_nod_icb_rigid_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & sph_bc_U%fdm2_fix_fld_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & ipol%i_velo, ipol%i_vort)
else
- call cal_sph_nod_icb_rigid_velo2(ipol%i_velo)
- call cal_sph_nod_icb_rigid_rot2(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_icb_rigid_velo2(nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%r_ICB, vt_ICB_bc, ipol%i_velo)
+ call cal_sph_nod_icb_rigid_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & sph_bc_U%fdm2_fix_fld_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & ipol%i_velo, ipol%i_vort)
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)
else
- call cal_sph_nod_cmb_rigid_v_and_w(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_cmb_rigid_v_and_w &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & sph_bc_U%fdm2_fix_fld_CMB, sph_bc_U%fdm2_fix_dr_CMB, &
+ & vt_CMB_bc, ipol%i_velo, ipol%i_vort)
end if
!
call cal_sph_diff_pol_and_rot2(nlayer_ICB, nlayer_CMB, &
@@ -235,15 +247,19 @@
if (sph_bc_U%iflag_icb .eq. iflag_free_slip) then
call cal_sph_nod_icb_free_vpol2(i_field)
else if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then
- call cal_sph_nod_icb_rotate_velo2(i_field)
+ call cal_sph_nod_icb_rotate_velo2 &
+ & (idx_rj_degree_zero, idx_rj_degree_one, nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%r_ICB, vt_ICB_bc, i_field)
else
- call cal_sph_nod_icb_rigid_velo2(i_field)
+ call cal_sph_nod_icb_rigid_velo2(nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%r_ICB, vt_ICB_bc, i_field)
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
call cal_sph_nod_cmb_free_vpol2(i_field)
else
- call cal_sph_nod_cmb_rigid_velo2(i_field)
+ 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)
end if
!
call cal_sph_diff_poloidal2(nlayer_ICB, nlayer_CMB, 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 310e637..39f77f0 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
@@ -77,15 +77,24 @@
if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then
call cal_sph_nod_icb_free_rot2(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(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_icb_rigid_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & sph_bc_U%fdm2_fix_fld_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & ipol%i_velo, ipol%i_vort)
else
- call cal_sph_nod_icb_rigid_rot2(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_icb_rigid_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & sph_bc_U%fdm2_fix_fld_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & ipol%i_velo, ipol%i_vort)
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)
else
- call cal_sph_nod_cmb_rigid_rot2(ipol%i_velo, ipol%i_vort)
+ call cal_sph_nod_cmb_rigid_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & sph_bc_U%fdm2_fix_fld_CMB, sph_bc_U%fdm2_fix_dr_CMB, &
+ & ipol%i_velo, ipol%i_vort)
end if
!
call cal_sph_nod_vect_rot2(nlayer_ICB, nlayer_CMB, &
@@ -223,13 +232,19 @@
if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then
call cal_sph_nod_icb_free_rot2(is_fld, is_rot)
else
- call cal_sph_nod_icb_rigid_rot2(is_fld, is_rot)
+ call cal_sph_nod_icb_rigid_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & sph_bc_U%fdm2_fix_fld_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & is_fld, is_rot)
end if
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
call cal_sph_nod_cmb_free_rot2(is_fld, is_rot)
else
- call cal_sph_nod_cmb_rigid_rot2(is_fld, is_rot)
+ call cal_sph_nod_cmb_rigid_rot2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & sph_bc_U%fdm2_fix_fld_CMB, sph_bc_U%fdm2_fix_dr_CMB, &
+ & is_fld, is_rot)
end if
!
call cal_sph_nod_vect_w_div_rot2(nlayer_ICB, nlayer_CMB, &
@@ -260,8 +275,10 @@
call cal_sph_nod_icb_free_diffuse2(coef_d_velo, &
& ipol%i_velo, ipol%i_v_diffuse)
else
- call cal_sph_nod_icb_rigid_diffuse2(coef_d_velo, &
- & ipol%i_velo, ipol%i_v_diffuse)
+ call cal_sph_nod_icb_rigid_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_in, sph_bc_U%r_ICB, &
+ & sph_bc_U%fdm2_fix_fld_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & coef_d_velo, ipol%i_velo, ipol%i_v_diffuse)
end if
call cal_dsdr_sph_no_bc_in_2(nidx_rj(2), &
& sph_bc_U%kr_in, sph_bc_U%fdm2_fix_fld_ICB, &
@@ -276,8 +293,10 @@
call cal_sph_nod_cmb_free_diffuse2(coef_d_velo, &
& ipol%i_velo, ipol%i_v_diffuse)
else
- call cal_sph_nod_cmb_rigid_diffuse2(coef_d_velo, ipol%i_velo, &
- & ipol%i_v_diffuse)
+ call cal_sph_nod_cmb_rigid_diffuse2 &
+ & (nidx_rj(2), sph_bc_U%kr_out, sph_bc_U%r_CMB, &
+ & sph_bc_U%fdm2_fix_fld_CMB, sph_bc_U%fdm2_fix_dr_CMB, &
+ & coef_d_velo, ipol%i_velo, ipol%i_v_diffuse)
end if
call cal_dsdr_sph_no_bc_out_2(nidx_rj(2), &
& sph_bc_U%kr_out, sph_bc_U%fdm2_fix_fld_CMB, &
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90
index 43b53f2..e360dc3 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90
@@ -7,13 +7,18 @@
!>@brief Evaluate velocity with non-slip boundary at CMB
!!
!!@verbatim
-!! subroutine cal_sph_nod_cmb_rigid_v_and_w(is_fld, is_rot)
-!! subroutine cal_sph_nod_cmb_rigid_velo2(is_fld)
-!! subroutine cal_sph_nod_cmb_rigid_rot2(is_fld, is_rot)
-!! subroutine cal_sph_nod_cmb_rigid_diffuse2(coef_d, &
-!! & is_fld, is_diffuse)
-!! subroutine cal_sph_nod_cmb_rgd_w_diffuse2(coef_d, &
+!! subroutine cal_sph_nod_cmb_rigid_v_and_w(jmax, kr_out, r_CMB, &
+!! & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, Vt_CMB, &
+!! & is_fld, is_rot)
+!! subroutine cal_sph_nod_cmb_rigid_velo2(jmax, kr_out, r_CMB, &
+!! & Vt_CMB, is_fld)
+!! subroutine cal_sph_nod_cmb_rigid_rot2(jmax, kr_out, r_CMB, &
+!! & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, is_fld, is_rot)
+!! subroutine cal_sph_nod_cmb_rigid_diffuse2(jmax, kr_out, r_CMB, &
+!! & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, coef_d, &
!! & is_fld, is_diffuse)
+!! subroutine cal_sph_nod_cmb_rgd_w_diffuse2(jmax, kr_out, r_CMB, &
+!! & fdm2_fix_fld_CMB, coef_d, is_fld, is_diffuse)
!!@endverbatim
!!
!!@n @param coef_d Coefficient for diffusion term
@@ -26,11 +31,8 @@
use m_precision
!
use m_constants
- use m_boundary_params_sph_MHD
- use m_spheric_parameter
use m_schmidt_poly_on_rtm
use m_sph_spectr_data
- use m_coef_fdm_fixed_CMB
!
implicit none
!
@@ -40,29 +42,35 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_rigid_v_and_w(is_fld, is_rot)
+ subroutine cal_sph_nod_cmb_rigid_v_and_w(jmax, kr_out, r_CMB, &
+ & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, Vt_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) :: r_CMB(0:2)
+ 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), intent(in) :: Vt_CMB(jmax)
!
integer(kind = kint) :: inod, j, i_n1, i_n2
real(kind = kreal) :: d2s_dr2, d1t_dr1
!
!
!$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
!
d_rj(inod,is_fld ) = zero
d_rj(inod,is_fld+1) = zero
- d_rj(inod,is_fld+2) = vt_CMB_bc(j)*ar_1d_rj(nlayer_CMB,1)
+ d_rj(inod,is_fld+2) = Vt_CMB(j)*r_CMB(1)
!
- d2s_dr2 = coef_fdm_fix_dr_CMB_2(-1,3) * d_rj(i_n1,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)
+ d2s_dr2 = fdm2_fix_dr_CMB(-1,3) * d_rj(i_n1,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
@@ -74,20 +82,24 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_rigid_velo2(is_fld)
+ subroutine cal_sph_nod_cmb_rigid_velo2(jmax, kr_out, r_CMB, &
+ & Vt_CMB, is_fld)
!
+ integer(kind = kint), intent(in) :: jmax, kr_out
integer(kind = kint), intent(in) :: is_fld
+ real(kind = kreal), intent(in) :: r_CMB(0:2)
+ real(kind = kreal), intent(in) :: Vt_CMB(jmax)
!
integer(kind = kint) :: inod, j
!
!
!$omp parallel do private(inod)
- do j = 1, nidx_rj(2)
- inod = j + (nlayer_CMB-1) * nidx_rj(2)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
!
d_rj(inod,is_fld ) = zero
d_rj(inod,is_fld+1) = zero
- d_rj(inod,is_fld+2) = vt_CMB_bc(j)*ar_1d_rj(nlayer_CMB,1)
+ d_rj(inod,is_fld+2) = Vt_CMB(j)*r_CMB(1)
end do
!$omp end parallel do
!
@@ -95,31 +107,35 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_rigid_rot2(is_fld, is_rot)
+ subroutine cal_sph_nod_cmb_rigid_rot2(jmax, kr_out, r_CMB, &
+ & 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) :: r_CMB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_fld_CMB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
!
integer(kind = kint) :: inod, j, i_n1, i_n2
real(kind = kreal) :: d2s_dr2, d1t_dr1
!
!
!$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_fix_dr_CMB_2(-1,3) * d_rj(i_n1,is_fld ) &
- & + coef_fdm_fix_dr_CMB_2( 0,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)
+ d2s_dr2 = fdm2_fix_dr_CMB(-1,3) * d_rj(i_n1,is_fld ) &
+ & + fdm2_fix_dr_CMB( 0,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 ) )
+ & * r_CMB(2)*d_rj(inod,is_fld ) )
end do
!$omp end parallel do
!
@@ -127,33 +143,37 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_rigid_diffuse2(coef_d, &
+ subroutine cal_sph_nod_cmb_rigid_diffuse2(jmax, kr_out, r_CMB, &
+ & fdm2_fix_fld_CMB, fdm2_fix_dr_CMB, 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) :: fdm2_fix_fld_CMB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,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_fix_dr_CMB_2(-1,3) * d_rj(i_n1,is_fld ) &
- & + coef_fdm_fix_dr_CMB_2( 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)
+ d2s_dr2 = fdm2_fix_dr_CMB(-1,3) * d_rj(i_n1,is_fld ) &
+ & + fdm2_fix_dr_CMB( 0,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)*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
!
@@ -161,36 +181,36 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_cmb_rgd_w_diffuse2(coef_d, &
- & is_fld, is_diffuse)
+ subroutine cal_sph_nod_cmb_rgd_w_diffuse2(jmax, kr_out, r_CMB, &
+ & fdm2_fix_fld_CMB, 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) :: coef_d
+ real(kind = kreal), intent(in) :: r_CMB(0:2)
+ 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)
-!
- d2s_dr2 = coef_fdm_fix_CMB_2(2,3) * d_rj(i_n2,is_fld ) &
- & + coef_fdm_fix_CMB_2(1,3) * d_rj(i_n1,is_fld ) &
- & + coef_fdm_fix_CMB_2(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)
+ do j = 1, jmax
+ inod = j + (kr_out-1) * jmax
+ i_n1 = inod - jmax
+ i_n2 = i_n1 - jmax
+!
+ d2s_dr2 = fdm2_fix_fld_CMB(2,3) * d_rj(i_n2,is_fld ) &
+ & + fdm2_fix_fld_CMB(1,3) * d_rj(i_n1,is_fld ) &
+ & + fdm2_fix_fld_CMB(0,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)*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_rigid_ICB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_ICB.f90
index 8a59900..ef1c356 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
@@ -7,13 +7,17 @@
!>@brief Evaluate velocity with non-slip boundary at ICB
!!
!!@verbatim
-!! subroutine cal_sph_nod_icb_rigid_velo2(is_fld)
-!! subroutine cal_sph_nod_icb_rotate_velo2(is_fld)
-!! subroutine cal_sph_nod_icb_rigid_rot2(is_fld, is_rot)
-!! subroutine cal_sph_nod_icb_rigid_diffuse2(coef_d, &
-!! & is_fld, is_diffuse)
-!! subroutine cal_sph_nod_icb_rgd_w_diffuse2(coef_d, &
+!! subroutine cal_sph_nod_icb_rigid_velo2(jmax, kr_in, r_ICB, &
+!! & Vt_ICB, is_fld)
+!! subroutine cal_sph_nod_icb_rotate_velo2(idx_rj_degree_zero, &
+!! & idx_rj_degree_one, jmax, kr_in, r_ICB, Vt_ICB, is_fld)
+!! subroutine cal_sph_nod_icb_rigid_rot2(jmax, kr_in, r_ICB, &
+!! & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, is_fld, is_rot)
+!! subroutine cal_sph_nod_icb_rigid_diffuse2(jmax, kr_in, r_ICB, &
+!! & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, coef_d, &
!! & is_fld, is_diffuse)
+!! subroutine cal_sph_nod_icb_rgd_w_diffuse2(jmax, kr_in, r_ICB, &
+!! & fdm2_fix_fld_ICB, coef_d, is_fld, is_diffuse)
!!@endverbatim
!!
!!@n @param coef_d Coefficient for diffusion term
@@ -26,11 +30,8 @@
use m_precision
!
use m_constants
- use m_spheric_parameter
use m_schmidt_poly_on_rtm
use m_sph_spectr_data
- use m_boundary_params_sph_MHD
- use m_coef_fdm_fixed_ICB
!
implicit none
!
@@ -41,20 +42,24 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_rigid_velo2(is_fld)
+ subroutine cal_sph_nod_icb_rigid_velo2(jmax, kr_in, r_ICB, &
+ & Vt_ICB, is_fld)
!
+ integer(kind = kint), intent(in) :: jmax, kr_in
integer(kind = kint), intent(in) :: is_fld
+ real(kind = kreal), intent(in) :: Vt_ICB(jmax)
+ real(kind = kreal), intent(in) :: r_ICB(0:2)
integer(kind = kint) :: inod, j, k
!
!
!$omp parallel do private(k,inod)
- do j = 1, nidx_rj(2)
- do k = 1, nlayer_ICB
- inod = j + (k-1) * nidx_rj(2)
+ do j = 1, jmax
+ do k = 1, kr_in
+ inod = j + (k-1) * jmax
!
d_rj(inod,is_fld ) = zero
d_rj(inod,is_fld+1) = zero
- d_rj(inod,is_fld+2) = vt_ICB_bc(j)*ar_1d_rj(k,1)
+ d_rj(inod,is_fld+2) = Vt_ICB(j)*r_ICB(1)
end do
end do
!$omp end parallel do
@@ -63,17 +68,23 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_rotate_velo2(is_fld)
+ subroutine cal_sph_nod_icb_rotate_velo2(idx_rj_degree_zero, &
+ & idx_rj_degree_one, jmax, kr_in, r_ICB, Vt_ICB, is_fld)
!
+ integer(kind = kint), intent(in) :: idx_rj_degree_zero
+ integer(kind = kint), intent(in) :: idx_rj_degree_one(-1:1)
+ integer(kind = kint), intent(in) :: jmax, kr_in
integer(kind = kint), intent(in) :: is_fld
+ real(kind = kreal), intent(in) :: Vt_ICB(jmax)
+ real(kind = kreal), intent(in) :: r_ICB(0:2)
!
integer(kind = kint) :: inod, iICB, j, k
!
!
!$omp parallel do private(k,inod)
- do j = 1, nidx_rj(2)
- do k = 1, nlayer_ICB
- inod = j + (k-1) * nidx_rj(2)
+ do j = 1, jmax
+ do k = 1, kr_in
+ inod = j + (k-1) * jmax
!
d_rj(inod,is_fld ) = zero
d_rj(inod,is_fld+1) = zero
@@ -82,30 +93,28 @@
!$omp end parallel do
!
if(idx_rj_degree_zero .gt. 0) then
- do k = 1, nlayer_ICB
- inod = idx_rj_degree_zero + (k-1) * nidx_rj(2)
- d_rj(inod,is_fld+2) &
- & = vt_ICB_bc(idx_rj_degree_zero)*ar_1d_rj(k,1)
+ do k = 1, kr_in
+ inod = idx_rj_degree_zero + (k-1) * jmax
+ d_rj(inod,is_fld+2) = Vt_ICB(idx_rj_degree_zero)*r_ICB(1)
end do
end if
!
!$omp parallel do private(k,inod)
- do j = idx_rj_degree_one(1)+1, nidx_rj(2)
- do k = 1, nlayer_ICB
- inod = j + (k-1) * nidx_rj(2)
- d_rj(inod,is_fld+2) = vt_ICB_bc(j)*ar_1d_rj(k,1)
+ do j = idx_rj_degree_one(1)+1, jmax
+ do k = 1, kr_in
+ inod = j + (k-1) * jmax
+ d_rj(inod,is_fld+2) = Vt_ICB(j)*r_ICB(1)
end do
end do
!$omp end parallel do
!
do j = -1, 1
if(idx_rj_degree_one(j) .gt. 0) then
- iICB = idx_rj_degree_one(j) + (nlayer_ICB-1) * nidx_rj(2)
- do k = 1, nlayer_ICB-1
- inod = idx_rj_degree_one(j) + (k-1) * nidx_rj(2)
- d_rj(inod,is_fld+2) &
- & = d_rj(iICB,is_fld+2)*ar_1d_rj(nlayer_ICB,2) &
- & * radius_1d_rj_r(k)*radius_1d_rj_r(k)
+ iICB = idx_rj_degree_one(j) + (kr_in-1) * jmax
+ do k = 1, kr_in-1
+ inod = idx_rj_degree_one(j) + (k-1) * jmax
+ d_rj(inod,is_fld+2) = d_rj(iICB,is_fld+2)*r_ICB(2) &
+ & * r_ICB(0)*r_ICB(0)
end do
end if
end do
@@ -114,31 +123,36 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_rigid_rot2(is_fld, is_rot)
+ subroutine cal_sph_nod_icb_rigid_rot2(jmax, kr_in, r_ICB, &
+ & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, is_fld, is_rot)
!
+ integer(kind = kint), intent(in) :: jmax, kr_in
integer(kind = kint), intent(in) :: is_fld
integer(kind = kint), intent(in) :: is_rot
+ 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) :: fdm2_fix_dr_ICB(-1:1,3)
!
integer(kind = kint) :: inod, j, i_p1, i_p2
real(kind = kreal) :: d2s_dr2, d1t_dr1
!
!
!$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_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)
+ d2s_dr2 = 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 ))
+ d_rj(inod,is_rot+2) = - (d2s_dr2 &
+ & - g_sph_rj(j,3) * r_ICB(2)*d_rj(inod,is_fld ))
end do
!$omp end parallel do
!
@@ -146,35 +160,37 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_rigid_diffuse2(coef_d, &
+ subroutine cal_sph_nod_icb_rigid_diffuse2(jmax, kr_in, r_ICB, &
+ & fdm2_fix_fld_ICB, fdm2_fix_dr_ICB, 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) :: fdm2_fix_fld_ICB(0:2,3)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1: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_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)
+ d2s_dr2 = 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)*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
!
@@ -182,36 +198,36 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_sph_nod_icb_rgd_w_diffuse2(coef_d, &
- & is_fld, is_diffuse)
+ subroutine cal_sph_nod_icb_rgd_w_diffuse2(jmax, kr_in, r_ICB, &
+ & fdm2_fix_fld_ICB, 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) :: fdm2_fix_fld_ICB(0:2,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)
-!
- d2s_dr2 = coef_fdm_fix_ICB_2( 0,3) * d_rj(inod,is_fld ) &
- & + coef_fdm_fix_ICB_2( 1,3) * d_rj(i_p1,is_fld ) &
- & + coef_fdm_fix_ICB_2( 2,3) * d_rj(i_p2,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 = fdm2_fix_fld_ICB( 0,3) * d_rj(inod,is_fld ) &
+ & + fdm2_fix_fld_ICB( 1,3) * d_rj(i_p1,is_fld ) &
+ & + fdm2_fix_fld_ICB( 2,3) * d_rj(i_p2,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)*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
!
More information about the CIG-COMMITS
mailing list