[cig-commits] [commit] Hiro_latest: Add new structures for basic boundary conditions for future replacement (29555ed)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Nov 18 16:21:42 PST 2013
Repository : ssh://geoshell/calypso
On branch : Hiro_latest
Link : https://github.com/geodynamics/calypso/compare/93e9f8f974c7a247c8f02e54ec18de063f86c8fb...3c548304673360ddedd7d68c8095b3fb74a2b9ce
>---------------------------------------------------------------
commit 29555edd1a7044c496eff774502555cc35ba8eaa
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Tue Nov 12 09:04:39 2013 -0800
Add new structures for basic boundary conditions for future replacement
>---------------------------------------------------------------
29555edd1a7044c496eff774502555cc35ba8eaa
.../MHD_src/sph_MHD/cal_fdm_coefs_4_boundaries.f90 | 288 +++++++++++++++++++++
.../MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90 | 13 +-
.../MHD_src/sph_MHD/m_boundary_params_sph_MHD.f90 | 66 +++++
.../MHD_src/sph_MHD/m_coef_fdm_fixed_CMB.f90 | 82 ------
.../MHD_src/sph_MHD/m_coef_fdm_fixed_ICB.f90 | 81 ------
.../MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90 | 217 ++++++++++++++++
6 files changed, 580 insertions(+), 167 deletions(-)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_fdm_coefs_4_boundaries.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_fdm_coefs_4_boundaries.f90
new file mode 100644
index 0000000..e8700cf
--- /dev/null
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_fdm_coefs_4_boundaries.f90
@@ -0,0 +1,288 @@
+!>@file cal_fdm_coefs_4_boundaries.f90
+!!@brief module cal_fdm_coefs_4_boundaries
+!!
+!!@author H. Matsui
+!!@date Programmed in Jan., 2010
+!!@n Modified in Nov., 2013
+!
+!>@brief Obtain FDM matrix for basic boundary conditions at boundaries
+!!
+!!@verbatim
+!! subroutine cal_fdm2_coef_fix_fld_ICB(r_from_ICB, &
+!! & fdm2_fix_fld_ICB)
+!! subroutine cal_fdm2_coef_fix_df_ICB(r_from_ICB, &
+!! & fdm2_fix_dr_ICB)
+!!
+!! subroutine cal_fdm2_coef_fix_fld_CMB(r_from_CMB2, &
+!! & fdm2_fix_fld_CMB)
+!! subroutine cal_fdm2_coef_fix_df_CMB(r_from_CMB1, &
+!! & fdm2_fix_dr_CMB)
+!!
+!! Matrix for derivatives with fixed field
+!! at inner boundary of the shell
+!! dfdr = fdm2_fix_fld_ICB( 0,2) * d_rj(ICB )
+!! + fdm2_fix_fld_ICB( 1,2) * d_rj(ICB+1)
+!! + fdm2_fix_fld_ICB( 2,2) * d_rj(ICB+2)
+!! d2fdr2 = fdm2_fix_fld_ICB( 0,3) * d_rj(ICB )
+!! + fdm2_fix_fld_ICB( 1,3) * d_rj(ICB+1)
+!! + fdm2_fix_fld_ICB( 2,3) * d_rj(ICB+2)
+!!
+!! Matrix for field and 2nd derivatives with fixed gradient
+!! at inner boundary of the shell
+!! d_rj(k) = fdm2_fix_dr_ICB(-1,2) * dfdr(ICB)
+!! + fdm2_fix_dr_ICB( 0,2) * d_rj(ICB )
+!! + fdm2_fix_dr_ICB( 1,2) * d_rj(ICB+1)
+!! d2fdr2 = fdm2_fix_dr_ICB(-1,3) * dfdr(ICB)
+!! + fdm2_fix_dr_ICB( 0,3) * d_rj(ICB )
+!! + fdm2_fix_dr_ICB( 1,3) * d_rj(ICB+1)
+!!
+!! Matrix for derivatives with fixed field
+!! at outer boundary of the shell
+!! dfdr = fdm2_fix_fld_CMB( 2,2) * d_rj(CMB-2)
+!! + fdm2_fix_fld_CMB( 1,2) * d_rj(CMB-1)
+!! + fdm2_fix_fld_CMB( 0,2) * d_rj(CMB )
+!! d2fdr2 = fdm2_fix_fld_CMB( 2,3) * d_rj(CMB-2)
+!! + fdm2_fix_fld_CMB( 1,3) * d_rj(CMB-1)
+!! + fdm2_fix_fld_CMB( 0,3) * d_rj(CMB )
+!!
+!! Matrix for field and 2nd derivatives with fixed gradient
+!! at outer boundary of the shell
+!! d_rj(k) = fdm2_fix_dr_CMB(-1,1) * d_rj(CMB-1)
+!! + fdm2_fix_dr_CMB( 0,1) * d_rj(CMB )
+!! + fdm2_fix_dr_CMB( 1,1) * dfdr(CMB)
+!! d2fdr2 = fdm2_fix_dr_CMB(-1,3) * d_rj(CMB-1)
+!! + fdm2_fix_dr_CMB( 0,3) * d_rj(CMB )
+!! + fdm2_fix_dr_CMB( 1,3) * dfdr(CMB)
+!!@endverbatim
+!!
+!!@n @param r_from_ICB(0:2) radius to teo next points of ICB
+!!@n @param r_from_CMB2(-2:0) radius from two next points to CMB
+!!@n @param r_from_CMB1(-1:0) radius from next points to CMB
+!!
+ module cal_fdm_coefs_4_boundaries
+!
+ use m_precision
+!
+ use m_constants
+! use m_spheric_parameter
+ use cal_inverse_small_matrix
+!
+ implicit none
+!
+!> Work matrix to evaluate fdm2_fix_fld_ICB(0:2,3)
+!!@verbatim
+!! dfdr = mat_fdm_ICB_fix_2(2,1) * d_rj(ICB )
+!! + mat_fdm_ICB_fix_2(2,2) * d_rj(ICB+1)
+!! + mat_fdm_ICB_fix_2(2,3) * d_rj(ICB+2)
+!! d2fdr2 = mat_fdm_ICB_fix_2(3,1) * d_rj(ICB )
+!! + mat_fdm_ICB_fix_2(3,2) * d_rj(ICB+1)
+!! + mat_fdm_ICB_fix_2(3,3) * d_rj(ICB+2)
+!!@endverbatim
+ real(kind = kreal) :: mat_fdm_ICB_fix_2(3,3)
+!
+!> Work matrix to evaluate fdm2_fix_dr_ICB(-1:1,3)
+!!@verbatim
+!! d_rj(k) = mat_fdm_ICB_fix_dr_2(2,1) * d_rj(ICB )
+!! + mat_fdm_ICB_fix_dr_2(2,2) * dfdr(ICB)
+!! + mat_fdm_ICB_fix_dr_2(2,3) * d_rj(ICB+1)
+!! d2fdr2 = mat_fdm_ICB_fix_dr_2(3,1) * d_rj(ICB )
+!! + mat_fdm_ICB_fix_dr_2(3,2) * dfdr(ICB)
+!! + mat_fdm_ICB_fix_dr_2(3,3) * d_rj(ICB+1)
+!!@endverbatim
+ real(kind = kreal) :: mat_fdm_ICB_fix_dr_2(3,3)
+!
+!> Work matrix to evaluate fdm2_fix_fld_CMB(0:2,3)
+!!@verbatim
+!! dfdr = mat_fdm_CMB_fix_2(2,1) * d_rj(CMB )
+!! + mat_fdm_CMB_fix_2(2,2) * d_rj(CMB-1)
+!! + mat_fdm_CMB_fix_2(2,3) * d_rj(CMB-2)
+!! d2fdr2 = mat_fdm_CMB_fix_2(3,1) * d_rj(CMB )
+!! + mat_fdm_CMB_fix_2(3,2) * d_rj(CMB-1)
+!! + mat_fdm_CMB_fix_2(3,3) * d_rj(CMB-2)
+!!@endverbatim
+ real(kind = kreal) :: mat_fdm_CMB_fix_2(3,3)
+!
+!> Work matrix to evaluate fdm2_fix_dr_CMB(-1:1,3)
+!!@verbatim
+!! d_rj(k) = mat_fdm_CMB_fix_dr_2(1,1) * d_rj(CMB )
+!! + mat_fdm_CMB_fix_dr_2(1,2) * dfdr(CMB)
+!! + mat_fdm_CMB_fix_dr_2(1,3) * d_rj(CMB-1)
+!! d2fdr2 = mat_fdm_CMB_fix_dr_2(3,1) * d_rj(CMB )
+!! + mat_fdm_CMB_fix_dr_2(3,2) * dfdr(CMB)
+!! + mat_fdm_CMB_fix_dr_2(3,3) * d_rj(CMB-1)
+!!@endverbatim
+ real(kind = kreal) :: mat_fdm_CMB_fix_dr_2(3,3)
+!
+ private :: mat_fdm_ICB_fix_2, mat_fdm_ICB_fix_dr_2
+ private :: mat_fdm_CMB_fix_2, mat_fdm_CMB_fix_dr_2
+!
+! -----------------------------------------------------------------------
+!
+ contains
+!
+! -----------------------------------------------------------------------
+!
+ subroutine cal_fdm2_coef_fix_fld_ICB(r_from_ICB, &
+ & fdm2_fix_fld_ICB)
+!
+ real(kind = kreal), intent(in) :: r_from_ICB(0:2)
+ real(kind = kreal), intent(inout) :: fdm2_fix_fld_ICB(0:2,3)
+!
+ integer(kind = kint) :: ierr
+ real(kind = kreal) :: mat_taylor_3(3,3)
+ real(kind = kreal) :: dr_p1, dr_p2
+!
+!
+ dr_p1 = r_from_ICB(1) - r_from_ICB(0)
+ dr_p2 = r_from_ICB(2) - r_from_ICB(0)
+!
+ mat_taylor_3(1,1) = one
+ mat_taylor_3(1,2) = zero
+ mat_taylor_3(1,3) = zero
+!
+ mat_taylor_3(2,1) = one
+ mat_taylor_3(2,2) = dr_p1
+ mat_taylor_3(2,3) = dr_p1*dr_p1 / two
+!
+ mat_taylor_3(3,1) = one
+ mat_taylor_3(3,2) = dr_p2
+ mat_taylor_3(3,3) = dr_p2*dr_p2 / two
+!
+ call cal_inverse_33_matrix(mat_taylor_3, mat_fdm_ICB_fix_2, ierr)
+ if(ierr .eq. 1) then
+ write(*,*) 'singular matrix cal_fdm2_coef_fix_fld_ICB ', &
+ & r_from_ICB(0)
+ end if
+!
+ fdm2_fix_fld_ICB(0,1:3) = mat_fdm_ICB_fix_2(1:3,1)
+ fdm2_fix_fld_ICB(1,1:3) = mat_fdm_ICB_fix_2(1:3,2)
+ fdm2_fix_fld_ICB(2,1:3) = mat_fdm_ICB_fix_2(1:3,3)
+!
+ end subroutine cal_fdm2_coef_fix_fld_ICB
+!
+! -----------------------------------------------------------------------
+!
+ subroutine cal_fdm2_coef_fix_df_ICB(r_from_ICB, &
+ & fdm2_fix_dr_ICB)
+!
+ real(kind = kreal), intent(in) :: r_from_ICB(0:1)
+ real(kind = kreal), intent(inout) :: fdm2_fix_dr_ICB(-1:1,3)
+!
+ integer(kind = kint) :: ierr
+ real(kind = kreal) :: mat_taylor_3(3,3)
+ real(kind = kreal) :: dr_p1
+!
+!
+ dr_p1 = r_from_ICB(1) - r_from_ICB(0)
+!
+ mat_taylor_3(1,1) = one
+ mat_taylor_3(1,2) = zero
+ mat_taylor_3(1,3) = zero
+!
+ mat_taylor_3(2,1) = zero
+ mat_taylor_3(2,2) = one
+ mat_taylor_3(2,3) = zero
+!
+ mat_taylor_3(3,1) = dr_p1
+ mat_taylor_3(3,2) = one
+ mat_taylor_3(3,3) = dr_p1*dr_p1 / two
+!
+ call cal_inverse_33_matrix(mat_taylor_3, mat_fdm_ICB_fix_dr_2, &
+ & ierr)
+ if(ierr .eq. 1) then
+ write(*,*) 'singular matrix cal_fdm2_coef_fix_df_ICB ', &
+ & r_from_ICB(0)
+ end if
+!
+ fdm2_fix_dr_ICB(-1,1:3) = mat_fdm_ICB_fix_dr_2(1:3,1)
+ fdm2_fix_dr_ICB( 0,1:3) = mat_fdm_ICB_fix_dr_2(1:3,2)
+ fdm2_fix_dr_ICB( 1,1:3) = mat_fdm_ICB_fix_dr_2(1:3,3)
+!
+ end subroutine cal_fdm2_coef_fix_df_ICB
+!
+! -----------------------------------------------------------------------
+! -----------------------------------------------------------------------
+!
+ subroutine cal_fdm2_coef_fix_fld_CMB(r_from_CMB2, &
+ & fdm2_fix_fld_CMB)
+!
+ real(kind = kreal), intent(in) :: r_from_CMB2(-2:0)
+ real(kind = kreal), intent(inout) :: fdm2_fix_fld_CMB(0:2,3)
+!
+ integer(kind = kint) :: ierr
+ real(kind = kreal) :: mat_taylor_3(3,3)
+ real(kind = kreal) :: dr_n1, dr_n2
+!
+!
+ dr_n1 = r_from_CMB2(0) - r_from_CMB2(-1)
+ dr_n2 = r_from_CMB2(0) - r_from_CMB2(-2)
+!
+ mat_taylor_3(1,1) = one
+ mat_taylor_3(1,2) = zero
+ mat_taylor_3(1,3) = zero
+!
+ mat_taylor_3(2,1) = one
+ mat_taylor_3(2,2) =-dr_n1
+ mat_taylor_3(2,3) = dr_n1*dr_n1 / two
+!
+ mat_taylor_3(3,1) = one
+ mat_taylor_3(3,2) =-dr_n2
+ mat_taylor_3(3,3) = dr_n2*dr_n2 / two
+!
+ call cal_inverse_33_matrix(mat_taylor_3, mat_fdm_CMB_fix_2, ierr)
+!
+ if(ierr .eq. 1) then
+ write(*,*) 'singular matrix cal_fdm2_coef_fix_fld_CMB ', &
+ & r_from_CMB2(0)
+ end if
+!
+ fdm2_fix_fld_CMB(2,1:3) = mat_fdm_CMB_fix_2(1:3,3)
+ fdm2_fix_fld_CMB(1,1:3) = mat_fdm_CMB_fix_2(1:3,2)
+ fdm2_fix_fld_CMB(0,1:3) = mat_fdm_CMB_fix_2(1:3,1)
+!
+ end subroutine cal_fdm2_coef_fix_fld_CMB
+!
+! -----------------------------------------------------------------------
+!
+ subroutine cal_fdm2_coef_fix_df_CMB(r_from_CMB1, &
+ & fdm2_fix_dr_CMB)
+!
+ real(kind = kreal), intent(in) :: r_from_CMB1(-1:0)
+ real(kind = kreal), intent(inout) :: fdm2_fix_dr_CMB(-1:1,3)
+!
+ integer(kind = kint) :: ierr
+ real(kind = kreal) :: mat_taylor_3(3,3)
+ real(kind = kreal) :: dr_n1
+!
+!
+ dr_n1 = r_from_CMB1(0) - r_from_CMB1(-1)
+!
+ mat_taylor_3(1,1) = one
+ mat_taylor_3(1,2) = zero
+ mat_taylor_3(1,3) = zero
+!
+ mat_taylor_3(2,1) = zero
+ mat_taylor_3(2,2) = one
+ mat_taylor_3(2,3) = zero
+!
+ mat_taylor_3(3,1) = one
+ mat_taylor_3(3,2) =-dr_n1
+ mat_taylor_3(3,3) = dr_n1*dr_n1 / two
+!
+ call cal_inverse_33_matrix(mat_taylor_3, mat_fdm_CMB_fix_dr_2, &
+ & ierr)
+!
+ if(ierr .eq. 1) then
+ write(*,*) 'singular matrix cal_fdm2_coef_fix_df_CMB ', &
+ & r_from_CMB1(0)
+ end if
+!
+ fdm2_fix_dr_CMB(-1,1:3) = mat_fdm_CMB_fix_dr_2(1:3,3)
+ fdm2_fix_dr_CMB( 0,1:3) = mat_fdm_CMB_fix_dr_2(1:3,1)
+ fdm2_fix_dr_CMB( 1,1:3) = mat_fdm_CMB_fix_dr_2(1:3,2)
+!
+ end subroutine cal_fdm2_coef_fix_df_CMB
+!
+! -----------------------------------------------------------------------
+!
+ end module cal_fdm_coefs_4_boundaries
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90
index f623251..cf55107 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_bc_fdm_matrix.f90
@@ -32,13 +32,18 @@
use m_coef_fdm_free_ICB
use m_coef_fdm_free_CMB
use m_coef_fdm_to_center
+ use cal_fdm_coefs_4_boundaries
!
!
- call cal_2nd_nod_ICB_fixed_fdm(radius_1d_rj_r(nlayer_ICB))
- call cal_2nd_nod_ICB_fix_df_fdm(radius_1d_rj_r(nlayer_ICB))
+ call cal_fdm2_coef_fix_fld_ICB(radius_1d_rj_r(nlayer_ICB), &
+ & coef_fdm_fix_ICB_2)
+ call cal_fdm2_coef_fix_df_ICB(radius_1d_rj_r(nlayer_ICB), &
+ & coef_fdm_fix_dr_ICB_2)
!
- call cal_2nd_nod_CMB_fixed_fdm(radius_1d_rj_r(nlayer_CMB-2))
- call cal_2nd_nod_CMB_fix_df_fdm(radius_1d_rj_r(nlayer_CMB-1))
+ call cal_fdm2_coef_fix_fld_CMB(radius_1d_rj_r(nlayer_CMB-2), &
+ & coef_fdm_fix_CMB_2)
+ call cal_fdm2_coef_fix_df_CMB(radius_1d_rj_r(nlayer_CMB-1), &
+ & coef_fdm_fix_dr_CMB_2)
!
!
call cal_2nd_ICB_free_vp_bc_fdm(radius_1d_rj_r(nlayer_ICB))
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/m_boundary_params_sph_MHD.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/m_boundary_params_sph_MHD.f90
new file mode 100644
index 0000000..565dfad
--- /dev/null
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/m_boundary_params_sph_MHD.f90
@@ -0,0 +1,66 @@
+!>@file m_boundary_params_sph_MHD.f90
+!!@brief module m_boundary_params_sph_MHD
+!!
+!!@author H. Matsui
+!!@date Programmed in Oct., 2009
+!
+!>@brief Structure for basic boundary conditions for spherical dynamo
+!!
+!!
+!!@verbatim
+!! subroutine alloc_fixed_bc_array(jmax, bc_param)
+!! subroutine dealloc_fixed_bc_array(bc_param)
+!!
+!! subroutine cal_fdm_coefs_4_BCs(nri, radius, bc_param)
+!! subroutine check_fdm_coefs_4_BC2(label, bc_param)
+!! type(sph_boundary_type), intent(inout) :: bc_param
+!!@endverbatim
+!!
+!!@n @param jmax number of modes for spherical harmonics @f$L*(L+2)@f$
+!!@n @param nri number of radial grid points
+!!@n @param radius radius
+!
+ module m_boundary_params_sph_MHD
+!
+ use m_precision
+ use t_boundary_params_sph_MHD
+!
+ implicit none
+!
+!
+!> Structure for basic velocity boundary condition parameters
+ type(sph_boundary_type), save :: sph_bc_U
+!> Structure for basic magnetic boundary condition parameters
+ type(sph_boundary_type), save :: sph_bc_B
+!> Structure for basic thermal boundary condition parameters
+ type(sph_boundary_type), save :: sph_bc_T
+!> Structure for basic compositional boundary condition parameters
+ type(sph_boundary_type), save :: sph_bc_C
+!
+! -----------------------------------------------------------------------
+!
+ contains
+!
+! -----------------------------------------------------------------------
+!
+ subroutine set_radial_range_by_BC(iflag_icb_bc, kr_in, kr_out)
+!
+ use m_spheric_parameter
+!
+ integer(kind = kint), intent(in) :: iflag_icb_bc
+ integer(kind = kint), intent(inout) :: kr_in, kr_out
+!
+!
+ if (iflag_icb_bc .eq. iflag_sph_fill_center &
+ & .or. iflag_icb_bc .eq. iflag_sph_fix_center) then
+ kr_in = nlayer_2_center
+ else
+ kr_in = nlayer_ICB
+ end if
+ kr_out = nlayer_CMB
+!
+ end subroutine set_radial_range_by_BC
+!
+! -----------------------------------------------------------------------
+!
+ end module m_boundary_params_sph_MHD
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_fixed_CMB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_fixed_CMB.f90
index 6555f66..80d2bb0 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_fixed_CMB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_fixed_CMB.f90
@@ -7,9 +7,6 @@
!>@brief Matrix to evaluate radial derivative at CMB
!!
!!@verbatim
-!! subroutine cal_2nd_nod_CMB_fixed_fdm(r_from_CMB2)
-!! subroutine cal_2nd_nod_CMB_fix_df_fdm(r_from_CMB1)
-!!
!! subroutine check_coef_fdm_fix_dr_CMB
!!
!! Matrix for derivatives with fixed field
@@ -79,85 +76,6 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_2nd_nod_CMB_fixed_fdm(r_from_CMB2)
-!
- real(kind = kreal) :: r_from_CMB2(-2:0)
-!
- integer(kind = kint) :: ierr
- real(kind = kreal) :: mat_taylor_3(3,3)
- real(kind = kreal) :: dr_n1, dr_n2
-!
-!
- dr_n1 = r_from_CMB2(0) - r_from_CMB2(-1)
- dr_n2 = r_from_CMB2(0) - r_from_CMB2(-2)
-!
- mat_taylor_3(1,1) = one
- mat_taylor_3(1,2) = zero
- mat_taylor_3(1,3) = zero
-!
- mat_taylor_3(2,1) = one
- mat_taylor_3(2,2) =-dr_n1
- mat_taylor_3(2,3) = dr_n1*dr_n1 / two
-!
- mat_taylor_3(3,1) = one
- mat_taylor_3(3,2) =-dr_n2
- mat_taylor_3(3,3) = dr_n2*dr_n2 / two
-!
- call cal_inverse_33_matrix(mat_taylor_3, mat_fdm_CMB_fix_2, ierr)
-!
- if(ierr .eq. 1) then
- write(*,*) 'singular matrix cal_2nd_nod_CMB_fixed_fdm ', &
- & r_from_CMB2(0)
- end if
-!
- coef_fdm_fix_CMB_2(2,1:3) = mat_fdm_CMB_fix_2(1:3,3)
- coef_fdm_fix_CMB_2(1,1:3) = mat_fdm_CMB_fix_2(1:3,2)
- coef_fdm_fix_CMB_2(0,1:3) = mat_fdm_CMB_fix_2(1:3,1)
-!
- end subroutine cal_2nd_nod_CMB_fixed_fdm
-!
-! -----------------------------------------------------------------------
-!
- subroutine cal_2nd_nod_CMB_fix_df_fdm(r_from_CMB1)
-!
- real(kind = kreal) :: r_from_CMB1(-1:0)
-!
- integer(kind = kint) :: ierr
- real(kind = kreal) :: mat_taylor_3(3,3)
- real(kind = kreal) :: dr_n1
-!
-!
- dr_n1 = r_from_CMB1(0) - r_from_CMB1(-1)
-!
- mat_taylor_3(1,1) = one
- mat_taylor_3(1,2) = zero
- mat_taylor_3(1,3) = zero
-!
- mat_taylor_3(2,1) = zero
- mat_taylor_3(2,2) = one
- mat_taylor_3(2,3) = zero
-!
- mat_taylor_3(3,1) = one
- mat_taylor_3(3,2) =-dr_n1
- mat_taylor_3(3,3) = dr_n1*dr_n1 / two
-!
- call cal_inverse_33_matrix(mat_taylor_3, mat_fdm_CMB_fix_dr_2, &
- & ierr)
-!
- if(ierr .eq. 1) then
- write(*,*) 'singular matrix cal_2nd_nod_CMB_fix_df_fdm ', &
- & r_from_CMB1(0)
- end if
-!
- coef_fdm_fix_dr_CMB_2(-1,1:3) = mat_fdm_CMB_fix_dr_2(1:3,3)
- coef_fdm_fix_dr_CMB_2( 0,1:3) = mat_fdm_CMB_fix_dr_2(1:3,1)
- coef_fdm_fix_dr_CMB_2( 1,1:3) = mat_fdm_CMB_fix_dr_2(1:3,2)
-!
- end subroutine cal_2nd_nod_CMB_fix_df_fdm
-!
-! -----------------------------------------------------------------------
-! -----------------------------------------------------------------------
-!
subroutine check_coef_fdm_fix_dr_CMB
!
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_fixed_ICB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_fixed_ICB.f90
index bfab75f..c9bec4e 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_fixed_ICB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/m_coef_fdm_fixed_ICB.f90
@@ -7,9 +7,6 @@
!>@brief Matrix to evaluate radial derivative at ICB
!!
!!@verbatim
-!! subroutine cal_2nd_nod_ICB_fixed_fdm(r_from_ICB)
-!! subroutine cal_2nd_nod_ICB_fix_df_fdm(r_from_ICB)
-!!
!! subroutine check_coef_fdm_fix_dr_ICB
!!
!! Matrix for derivatives with fixed field
@@ -36,7 +33,6 @@
use m_precision
!
use m_constants
-! use m_spheric_parameter
use cal_inverse_small_matrix
!
implicit none
@@ -78,83 +74,6 @@
!
! -----------------------------------------------------------------------
!
- subroutine cal_2nd_nod_ICB_fixed_fdm(r_from_ICB)
-!
- real(kind = kreal) :: r_from_ICB(0:2)
-!
- integer(kind = kint) :: ierr
- real(kind = kreal) :: mat_taylor_3(3,3)
- real(kind = kreal) :: dr_p1, dr_p2
-!
-!
- dr_p1 = r_from_ICB(1) - r_from_ICB(0)
- dr_p2 = r_from_ICB(2) - r_from_ICB(0)
-!
- mat_taylor_3(1,1) = one
- mat_taylor_3(1,2) = zero
- mat_taylor_3(1,3) = zero
-!
- mat_taylor_3(2,1) = one
- mat_taylor_3(2,2) = dr_p1
- mat_taylor_3(2,3) = dr_p1*dr_p1 / two
-!
- mat_taylor_3(3,1) = one
- mat_taylor_3(3,2) = dr_p2
- mat_taylor_3(3,3) = dr_p2*dr_p2 / two
-!
- call cal_inverse_33_matrix(mat_taylor_3, mat_fdm_ICB_fix_2, ierr)
- if(ierr .eq. 1) then
- write(*,*) 'singular matrix cal_2nd_nod_ICB_fixed_fdm ', &
- & r_from_ICB(0)
- end if
-!
- coef_fdm_fix_ICB_2(0,1:3) = mat_fdm_ICB_fix_2(1:3,1)
- coef_fdm_fix_ICB_2(1,1:3) = mat_fdm_ICB_fix_2(1:3,2)
- coef_fdm_fix_ICB_2(2,1:3) = mat_fdm_ICB_fix_2(1:3,3)
-!
- end subroutine cal_2nd_nod_ICB_fixed_fdm
-!
-! -----------------------------------------------------------------------
-!
- subroutine cal_2nd_nod_ICB_fix_df_fdm(r_from_ICB)
-!
- real(kind = kreal) :: r_from_ICB(0:1)
-!
- integer(kind = kint) :: ierr
- real(kind = kreal) :: mat_taylor_3(3,3)
- real(kind = kreal) :: dr_p1
-!
-!
- dr_p1 = r_from_ICB(1) - r_from_ICB(0)
-!
- mat_taylor_3(1,1) = one
- mat_taylor_3(1,2) = zero
- mat_taylor_3(1,3) = zero
-!
- mat_taylor_3(2,1) = zero
- mat_taylor_3(2,2) = one
- mat_taylor_3(2,3) = zero
-!
- mat_taylor_3(3,1) = dr_p1
- mat_taylor_3(3,2) = one
- mat_taylor_3(3,3) = dr_p1*dr_p1 / two
-!
- call cal_inverse_33_matrix(mat_taylor_3, mat_fdm_ICB_fix_dr_2, &
- & ierr)
- if(ierr .eq. 1) then
- write(*,*) 'singular matrix cal_2nd_nod_ICB_fix_df_fdm ', &
- & r_from_ICB(0)
- end if
-!
- coef_fdm_fix_dr_ICB_2(-1,1:3) = mat_fdm_ICB_fix_dr_2(1:3,1)
- coef_fdm_fix_dr_ICB_2( 0,1:3) = mat_fdm_ICB_fix_dr_2(1:3,2)
- coef_fdm_fix_dr_ICB_2( 1,1:3) = mat_fdm_ICB_fix_dr_2(1:3,3)
-!
- end subroutine cal_2nd_nod_ICB_fix_df_fdm
-!
-! -----------------------------------------------------------------------
-! -----------------------------------------------------------------------
-!
subroutine check_coef_fdm_fix_dr_ICB
!
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90
new file mode 100644
index 0000000..15e3f5f
--- /dev/null
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/t_boundary_params_sph_MHD.f90
@@ -0,0 +1,217 @@
+!>@file t_boundary_params_sph_MHD.f90
+!!@brief module t_boundary_params_sph_MHD
+!!
+!!@author H. Matsui
+!!@date Programmed in Oct., 2009
+!
+!>@brief Structure for basic boundary conditions for spherical dynamo
+!!
+!!
+!!@verbatim
+!! subroutine alloc_fixed_bc_array(jmax, sph_bc)
+!! subroutine dealloc_fixed_bc_array(sph_bc)
+!!
+!! subroutine cal_fdm_coefs_4_BCs(nri, radius, sph_bc)
+!! subroutine check_fdm_coefs_4_BC2(label, sph_bc)
+!! subroutine check_sph_boundary_spectra(label, jmax, j_rj, sph_bc)
+!! type(sph_boundary_type), intent(inout) :: sph_bc
+!!@endverbatim
+!!
+!!@n @param jmax number of modes for spherical harmonics @f$L*(L+2)@f$
+!!@n @param nri number of radial grid points
+!!@n @param radius radius
+!!@n @param j_rj local spherical harmionics modes for f(r,j)
+!
+ module t_boundary_params_sph_MHD
+!
+ use m_precision
+!
+ implicit none
+!
+!
+!> integer flag for fixed velocity boundary
+ integer(kind = kint), parameter :: iflag_fixed_field = 0
+!> integer flag for free-slip boundary
+ integer(kind = kint), parameter :: iflag_fixed_flux = 1
+!> integer flag for whole sphere model
+ integer(kind = kint), parameter :: iflag_sph_fill_center = 41
+!> integer flag for whole sphere model
+ integer(kind = kint), parameter :: iflag_sph_fix_center = 42
+!
+!
+!> Structure for basic boundary condition parameters
+ type sph_boundary_type
+!> boundary condition flag at ICB
+ integer(kind = kint) :: iflag_icb = iflag_fixed_field
+!> boundary condition flag at CMB
+ integer(kind = kint) :: iflag_cmb = iflag_fixed_field
+!
+!> Start radial address of fluid shell for @f$ f(r,j) @f$
+ integer(kind = kint) :: kr_in = 1
+!> End radial address of fluid shell for @f$ f(r,j) @f$
+ integer(kind = kint) :: kr_out = 1
+!
+!> Fixed composition spectrum for ICB
+ real(kind= kreal), allocatable :: ICB_fld(:)
+!> Fixed composition flux spectrum for ICB
+ real(kind= kreal), allocatable :: ICB_flux(:)
+!> Fixed composition spectrum for CMB
+ real(kind= kreal), allocatable :: CMB_fld(:)
+!> Fixed composition flux spectrum for CMB
+ real(kind= kreal), allocatable :: CMB_flux(:)
+!
+!> Matrix to evaluate radial derivative at ICB with fiexed field
+ real(kind = kreal) :: fdm2_fix_fld_ICB(0:2,3)
+!> Matrix to evaluate field at ICB with fiexed radial derivative
+ real(kind = kreal) :: fdm2_fix_dr_ICB(-1:1,3)
+!> Matrix to evaluate radial derivative at CMB with fiexed field
+ real(kind = kreal) :: fdm2_fix_fld_CMB(0:2,3)
+!> Matrix to evaluate field at CMB with fiexed radial derivative
+ real(kind = kreal) :: fdm2_fix_dr_CMB(-1:1,3)
+ end type sph_boundary_type
+!
+! -----------------------------------------------------------------------
+!
+ contains
+!
+! -----------------------------------------------------------------------
+!
+ subroutine alloc_fixed_bc_array(jmax, sph_bc)
+!
+ integer(kind= kint), intent(in) :: jmax
+ type(sph_boundary_type), intent(inout) :: sph_bc
+!
+ allocate(sph_bc%ICB_fld(jmax))
+ allocate(sph_bc%CMB_fld(jmax))
+ allocate(sph_bc%ICB_flux(jmax))
+ allocate(sph_bc%CMB_flux(jmax))
+ sph_bc%ICB_fld = 0.0d0
+ sph_bc%CMB_fld = 0.0d0
+ sph_bc%ICB_flux = 0.0d0
+ sph_bc%CMB_flux = 0.0d0
+!
+ end subroutine alloc_fixed_bc_array
+!
+! -----------------------------------------------------------------------
+!
+ subroutine dealloc_fixed_bc_array(sph_bc)
+!
+ type(sph_boundary_type), intent(inout) :: sph_bc
+!
+!
+ deallocate(sph_bc%ICB_fld, sph_bc%CMB_fld)
+ deallocate(sph_bc%ICB_flux, sph_bc%CMB_flux)
+!
+ end subroutine dealloc_fixed_bc_array
+!
+! -----------------------------------------------------------------------
+! -----------------------------------------------------------------------
+!
+ subroutine cal_fdm_coefs_4_BCs(nri, radius, sph_bc)
+!
+ use cal_fdm_coefs_4_boundaries
+!
+ integer(kind = kint), intent(in) :: nri
+ real(kind = kreal), intent(in) :: radius(nri)
+ type(sph_boundary_type), intent(inout) :: sph_bc
+!
+!
+ call cal_fdm2_coef_fix_fld_ICB(radius(sph_bc%kr_in), &
+ & sph_bc%fdm2_fix_fld_ICB)
+ call cal_fdm2_coef_fix_df_ICB(radius(sph_bc%kr_in), &
+ & sph_bc%fdm2_fix_dr_ICB)
+!
+ call cal_fdm2_coef_fix_fld_CMB(radius(sph_bc%kr_out-2), &
+ & sph_bc%fdm2_fix_fld_ICB)
+ call cal_fdm2_coef_fix_df_CMB(radius(sph_bc%kr_out-1), &
+ & sph_bc%fdm2_fix_dr_ICB)
+!
+ end subroutine cal_fdm_coefs_4_BCs
+!
+! -----------------------------------------------------------------------
+!
+ subroutine check_fdm_coefs_4_BC2(label, sph_bc)
+!
+ character(len=kchara), intent(in) :: label
+ type(sph_boundary_type), intent(in) :: sph_bc
+!
+!
+ write(50,*) ' Boundary condition matrix for ', trim(label)
+!
+ write(50,*) ' fdm2_fix_fld_ICB'
+ write(50,*) ' mat_fdm21, mat_fdm22, mat_fdm23'
+ write(50,'(1p9E25.15e3)') sph_bc%fdm2_fix_fld_ICB(0:2,2)
+ write(50,*) ' mat_fdm31, mat_fdm32, mat_fdm33'
+ write(50,'(1p9E25.15e3)') sph_bc%fdm2_fix_fld_ICB(0:2,3)
+!
+ write(50,*) ' fdm2_fix_dr_ICB'
+ write(50,*) ' mat_fdm21, mat_fdm22, mat_fdm23'
+ write(50,'(1p9E25.15e3)') sph_bc%fdm2_fix_dr_ICB(-1:1,2)
+ write(50,*) ' mat_fdm31, mat_fdm32, mat_fdm33'
+ write(50,'(1p9E25.15e3)') sph_bc%fdm2_fix_dr_ICB(-1:1,3)
+!
+!
+ write(50,*) ' fdm2_fix_fld_CMB'
+ write(50,*) ' mat_fdm21, mat_fdm22, mat_fdm23'
+ write(50,'(1p9E25.15e3)') sph_bc%fdm2_fix_fld_CMB(0:2,2)
+ write(50,*) ' mat_fdm31, mat_fdm32, mat_fdm33'
+ write(50,'(1p9E25.15e3)') sph_bc%fdm2_fix_fld_CMB(0:2,3)
+!
+ write(50,*) ' fdm2_fix_dr_CMB'
+ write(50,*) ' mat_fdm11, mat_fdm12, mat_fdm13'
+ write(50,'(1p9E25.15e3)') sph_bc%fdm2_fix_dr_CMB(-1:1,1)
+ write(50,*) ' mat_fdm31, mat_fdm32, mat_fdm33'
+ write(50,'(1p9E25.15e3)') sph_bc%fdm2_fix_dr_CMB(-1:1,3)
+!
+ end subroutine check_fdm_coefs_4_BC2
+!
+! -----------------------------------------------------------------------
+!
+ subroutine check_sph_boundary_spectra(label, jmax, j_rj, sph_bc)
+!
+ use m_spheric_parameter
+ use m_bc_data_list
+ use m_surf_data_list
+!
+ integer(kind = kint), intent(in) :: jmax
+ integer(kind = kint), intent(in) :: j_rj(jmax,3)
+ character(len=kchara), intent(in) :: label
+ type(sph_boundary_type), intent(in) :: sph_bc
+!
+ integer(kind = kint) :: j
+!
+!
+ write(50,*) ' Boundary condition spectra for ', trim(label)
+!
+ write(50,*) 'iflag_icb', sph_bc%iflag_icb
+ if(sph_bc%iflag_icb .eq. iflag_fixed_field) then
+ write(50,*) 'field at ICB '
+ do j = 1, jmax
+ write(50,*) j_rj(j,1:3), sph_bc%ICB_fld(j)
+ end do
+ end if
+ if(sph_bc%iflag_icb .eq. iflag_fixed_flux) then
+ write(50,*) 'flux at ICB '
+ do j = 1, jmax
+ write(50,*) j_rj(j,1:3), sph_bc%ICB_flux(j)
+ end do
+ end if
+!
+ write(50,*) 'iflag_cmb', sph_bc%iflag_cmb
+ if(sph_bc%iflag_cmb .eq. iflag_fixed_field) then
+ do j = 1, jmax
+ write(50,*) j_rj(j,1:3), sph_bc%CMB_fld(j)
+ end do
+ end if
+ if(sph_bc%iflag_cmb .eq. iflag_fixed_flux) then
+ write(50,*) 'flux at CMB '
+ do j = 1, jmax
+ write(50,*) j_rj(j,1:3), sph_bc%CMB_flux(j)
+ end do
+ end if
+!
+ end subroutine check_sph_boundary_spectra
+!
+! -----------------------------------------------------------------------
+!
+ end module t_boundary_params_sph_MHD
More information about the CIG-COMMITS
mailing list