Skip to content

Commit 3f5fd14

Browse files
hiromatsuihiromatsui
committed
Add new module for boundary condition selector for explicit scalar diffusion
1 parent cf08b8f commit 3f5fd14

File tree

6 files changed

+321
-289
lines changed

6 files changed

+321
-289
lines changed

‎src/Fortran_libraries/MHD_src/radial_FDM/sel_sph_exp_scl_boundaries.f90‎

Lines changed: 253 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,253 @@
1+
!> @file sel_sph_exp_scl_boundaries.f90
2+
!! module sel_sph_exp_scl_boundaries
3+
!!
4+
!! @author H. Matsui
5+
!! @date Programmed in Oct. 2009
6+
!
7+
!> @brief Evaluate radial delivatives
8+
!!
9+
!!@verbatim
10+
!! subroutine sel_ICB_sph_scalar_diffusion &
11+
!! & (sph_rj, sph_bc, ICB_Sspec, fdm2_center, &
12+
!! & g_sph_rj, coef_diffuse, scl_rj, dfs_rj)
13+
!! subroutine sel_ICB_sph_scalar_val_diffuse &
14+
!! & (sph_rj, sph_bc, ICB_Sspec, fdm2_center, g_sph_rj, &
15+
!! & coef_diffuse, k_ratio, dk_dr, scl_rj, dfs_rj)
16+
!! type(sph_rj_grid), intent(in) :: sph_rj
17+
!! type(fdm2_center_mat), intent(in) :: fdm2_center
18+
!! type(sph_boundary_type), intent(in) :: sph_bc
19+
!! type(sph_scalar_BC_coef), intent(in) :: ICB_Sspec
20+
!! real(kind = kreal), intent(in) :: coef_diffuse
21+
!! real(kind = kreal), intent(in) :: k_ratio(0:sph_rj%nidx_rj(1))
22+
!! real(kind = kreal), intent(in) :: dk_dr(0:sph_rj%nidx_rj(1))
23+
!! real(kind=kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
24+
!! real(kind = kreal), intent(inout) :: scl_rj(sph_rj%nnod_rj)
25+
!! real(kind = kreal), intent(inout) :: dfs_rj(sph_rj%nnod_rj)
26+
!!
27+
!! subroutine sel_CMB_sph_scalar_diffusion &
28+
!! & (sph_rj, sph_bc, CMB_Sspec, g_sph_rj, coef_diffuse, &
29+
!! & is_fld, is_diffuse, n_point, ntot_phys_rj, d_rj)
30+
!! subroutine sel_CMB_sph_scalar_val_diffuse &
31+
!! & (sph_rj, sph_bc, CMB_Sspec, g_sph_rj, &
32+
!! & coef_diffuse, k_ratio, dk_dr, scl_rj, dfs_rj)
33+
!! type(sph_rj_grid), intent(in) :: sph_rj
34+
!! type(sph_boundary_type), intent(in) :: sph_bc
35+
!! type(sph_scalar_BC_coef), intent(in) :: CMB_Sspec
36+
!! real(kind = kreal), intent(in) :: coef_diffuse
37+
!! real(kind = kreal), intent(in) :: k_ratio, dk_dr
38+
!! real(kind=kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
39+
!! real(kind = kreal), intent(inout) :: scl_rj(sph_rj%nnod_rj)
40+
!! real(kind = kreal), intent(inout) :: dfs_rj(sph_rj%nnod_rj)
41+
!!@endverbatim
42+
!!
43+
!!@param sph_bc Structure for basic boundary condition parameters
44+
!!@param coef_diffuse Diffusion coefficient
45+
!!
46+
!!@param is_fld Spherical hermonics data address for input vector
47+
!!@param is_diffuse Input spectr diffusiton term address
48+
!
49+
module sel_sph_exp_scl_boundaries
50+
!
51+
use m_precision
52+
use m_constants
53+
!
54+
use t_spheric_rj_data
55+
use t_boundary_params_sph_MHD
56+
use t_boundary_sph_spectr
57+
use t_coef_fdm2_centre
58+
!
59+
implicit none
60+
!
61+
! -----------------------------------------------------------------------
62+
!
63+
contains
64+
!
65+
! -----------------------------------------------------------------------
66+
!
67+
subroutine sel_ICB_sph_scalar_diffusion &
68+
& (sph_rj, sph_bc, ICB_Sspec, fdm2_center, &
69+
& g_sph_rj, coef_diffuse, scl_rj, dfs_rj)
70+
!
71+
use sph_exp_fix_scl_diffuse_ICB
72+
use sph_exp_fix_flx_diffuse_ICB
73+
use sph_filled_center_diffuse2
74+
use sph_fixed_center_diffuse2
75+
!
76+
type(sph_rj_grid), intent(in) :: sph_rj
77+
type(fdm2_center_mat), intent(in) :: fdm2_center
78+
type(sph_boundary_type), intent(in) :: sph_bc
79+
type(sph_scalar_BC_coef), intent(in) :: ICB_Sspec
80+
!
81+
real(kind = kreal), intent(in) :: coef_diffuse
82+
real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
83+
!
84+
real(kind = kreal), intent(inout) :: scl_rj(sph_rj%nnod_rj)
85+
real(kind = kreal), intent(inout) :: dfs_rj(sph_rj%nnod_rj)
86+
!
87+
!
88+
if (sph_bc%iflag_icb .eq. iflag_sph_fill_center) then
89+
call sph_filled_ctr_diffuse_ctr2 &
90+
& (sph_rj%inod_rj_center, sph_rj%idx_rj_degree_zero, &
91+
& sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, sph_bc%r_ICB, &
92+
& fdm2_center%dmat_fix_fld, fdm2_center%dmat_fix_dr, &
93+
& coef_diffuse, scl_rj, dfs_rj)
94+
else if(sph_bc%iflag_icb .eq. iflag_sph_fix_center) then
95+
call sph_fixed_ctr_diffuse_ctr1 &
96+
& (sph_rj%inod_rj_center, sph_rj%idx_rj_degree_zero, &
97+
& sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
98+
& sph_bc%r_ICB, fdm2_center%dmat_fix_fld, ICB_Sspec%S_BC, &
99+
& coef_diffuse, scl_rj, dfs_rj)
100+
else if(sph_bc%iflag_icb .eq. iflag_fixed_flux &
101+
& .or. sph_bc%iflag_icb .eq. iflag_evolve_flux) then
102+
call sph_in_fix_flux_scl_diffuse2 &
103+
& (sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
104+
& sph_bc%kr_in, sph_bc%r_ICB, sph_bc%fdm2_fix_dr_ICB, &
105+
& ICB_Sspec%S_BC, coef_diffuse, scl_rj, dfs_rj)
106+
! else if(sph_bc%iflag_icb .eq. iflag_fixed_field &
107+
! & .or. sph_bc%iflag_icb .eq. iflag_evolve_field) then
108+
else
109+
call sph_in_fix_scalar_diffuse2 &
110+
& (sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
111+
& sph_bc%kr_in, sph_bc%r_ICB, sph_bc%fdm2_fix_fld_ICB, &
112+
& ICB_Sspec%S_BC, coef_diffuse, scl_rj, dfs_rj)
113+
end if
114+
!
115+
end subroutine sel_ICB_sph_scalar_diffusion
116+
!
117+
! -----------------------------------------------------------------------
118+
!
119+
subroutine sel_ICB_sph_scalar_val_diffuse &
120+
& (sph_rj, sph_bc, ICB_Sspec, fdm2_center, g_sph_rj, &
121+
& coef_diffuse, k_ratio, dk_dr, scl_rj, dfs_rj)
122+
!
123+
use sph_exp_fix_scl_diffuse_ICB
124+
use sph_exp_fix_flx_diffuse_ICB
125+
use sph_filled_center_diffuse2
126+
use sph_fixed_center_diffuse2
127+
!
128+
type(sph_rj_grid), intent(in) :: sph_rj
129+
type(fdm2_center_mat), intent(in) :: fdm2_center
130+
type(sph_boundary_type), intent(in) :: sph_bc
131+
type(sph_scalar_BC_coef), intent(in) :: ICB_Sspec
132+
!
133+
real(kind = kreal), intent(in) :: coef_diffuse
134+
real(kind = kreal), intent(in) :: k_ratio(0:sph_rj%nidx_rj(1))
135+
real(kind = kreal), intent(in) :: dk_dr(0:sph_rj%nidx_rj(1))
136+
real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
137+
!
138+
real(kind = kreal), intent(inout) :: scl_rj(sph_rj%nnod_rj)
139+
real(kind = kreal), intent(inout) :: dfs_rj(sph_rj%nnod_rj)
140+
!
141+
!
142+
if (sph_bc%iflag_icb .eq. iflag_sph_fill_center) then
143+
call sph_filled_ctr_val_diffuse_ctr2 &
144+
& (sph_rj%inod_rj_center, sph_rj%idx_rj_degree_zero, &
145+
& sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, sph_bc%r_ICB, &
146+
& fdm2_center%dmat_fix_fld, fdm2_center%dmat_fix_dr, &
147+
& coef_diffuse, k_ratio(0), dk_dr(0), scl_rj, dfs_rj)
148+
else if(sph_bc%iflag_icb .eq. iflag_sph_fix_center) then
149+
call sph_fixed_ctr_val_diffuse_ctr1 &
150+
& (sph_rj%inod_rj_center, sph_rj%idx_rj_degree_zero, &
151+
& sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
152+
& sph_bc%r_ICB, fdm2_center%dmat_fix_fld, ICB_Sspec%S_BC, &
153+
& coef_diffuse, k_ratio(1), dk_dr(1), scl_rj, dfs_rj)
154+
else if(sph_bc%iflag_icb .eq. iflag_fixed_flux &
155+
& .or. sph_bc%iflag_icb .eq. iflag_evolve_flux) then
156+
call sph_in_fix_flux_val_diffuse2 &
157+
& (sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, sph_bc%kr_in, &
158+
& sph_bc%r_ICB, sph_bc%fdm2_fix_dr_ICB, ICB_Sspec%S_BC, &
159+
& coef_diffuse, k_ratio(sph_bc%kr_in), dk_dr(sph_bc%kr_in), &
160+
& scl_rj, dfs_rj)
161+
! else if(sph_bc%iflag_icb .eq. iflag_fixed_field &
162+
! & .or. sph_bc%iflag_icb .eq. iflag_evolve_field) then
163+
else
164+
call sph_in_fix_scl_val_diffuse2 &
165+
& (sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, sph_bc%kr_in, &
166+
& sph_bc%r_ICB, sph_bc%fdm2_fix_fld_ICB, ICB_Sspec%S_BC, &
167+
& coef_diffuse, k_ratio(sph_bc%kr_in), dk_dr(sph_bc%kr_in), &
168+
& scl_rj, dfs_rj)
169+
end if
170+
!
171+
end subroutine sel_ICB_sph_scalar_val_diffuse
172+
!
173+
! -----------------------------------------------------------------------
174+
! -----------------------------------------------------------------------
175+
!
176+
subroutine sel_CMB_sph_scalar_diffusion &
177+
& (sph_rj, sph_bc, CMB_Sspec, g_sph_rj, coef_diffuse, &
178+
& scl_rj, dfs_rj)
179+
!
180+
use sph_exp_fix_scl_diffuse_CMB
181+
use sph_exp_fix_flx_diffuse_CMB
182+
!
183+
type(sph_rj_grid), intent(in) :: sph_rj
184+
type(sph_boundary_type), intent(in) :: sph_bc
185+
type(sph_scalar_BC_coef), intent(in) :: CMB_Sspec
186+
!
187+
real(kind = kreal), intent(in) :: coef_diffuse
188+
real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
189+
!
190+
real(kind = kreal), intent(inout) :: scl_rj(sph_rj%nnod_rj)
191+
real(kind = kreal), intent(inout) :: dfs_rj(sph_rj%nnod_rj)
192+
!
193+
!
194+
if (sph_bc%iflag_cmb .eq. iflag_fixed_flux &
195+
& .or. sph_bc%iflag_cmb .eq. iflag_evolve_flux) then
196+
call sph_out_fix_flux_scl_diffuse2 &
197+
& (sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
198+
& sph_bc%kr_out, sph_bc%r_CMB, sph_bc%fdm2_fix_dr_CMB, &
199+
& CMB_Sspec%S_BC, coef_diffuse, scl_rj, dfs_rj)
200+
! else if(sph_bc%iflag_cmb .eq. iflag_fixed_field &
201+
! & .or. sph_bc%iflag_cmb .eq. iflag_evolve_field) then
202+
else
203+
call sph_out_fix_scalar_diffuse2 &
204+
& (sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
205+
& sph_bc%kr_out, sph_bc%r_CMB, sph_bc%fdm2_fix_fld_CMB, &
206+
& CMB_Sspec%S_BC, coef_diffuse, scl_rj, dfs_rj)
207+
end if
208+
!
209+
end subroutine sel_CMB_sph_scalar_diffusion
210+
!
211+
! -----------------------------------------------------------------------
212+
!
213+
subroutine sel_CMB_sph_scalar_val_diffuse &
214+
& (sph_rj, sph_bc, CMB_Sspec, g_sph_rj, &
215+
& coef_diffuse, k_ratio, dk_dr, scl_rj, dfs_rj)
216+
!
217+
use sph_exp_fix_scl_diffuse_CMB
218+
use sph_exp_fix_flx_diffuse_CMB
219+
!
220+
type(sph_rj_grid), intent(in) :: sph_rj
221+
type(sph_boundary_type), intent(in) :: sph_bc
222+
type(sph_scalar_BC_coef), intent(in) :: CMB_Sspec
223+
!
224+
real(kind = kreal), intent(in) :: coef_diffuse
225+
real(kind = kreal), intent(in) :: k_ratio, dk_dr
226+
real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
227+
!
228+
real(kind = kreal), intent(inout) :: scl_rj(sph_rj%nnod_rj)
229+
real(kind = kreal), intent(inout) :: dfs_rj(sph_rj%nnod_rj)
230+
!
231+
!
232+
if (sph_bc%iflag_cmb .eq. iflag_fixed_flux &
233+
& .or. sph_bc%iflag_cmb .eq. iflag_evolve_flux) then
234+
call sph_out_fix_flux_val_diffuse2 &
235+
& (sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
236+
& sph_bc%kr_out, sph_bc%r_CMB, sph_bc%fdm2_fix_dr_CMB, &
237+
& CMB_Sspec%S_BC, coef_diffuse, k_ratio, dk_dr, &
238+
& scl_rj, dfs_rj)
239+
! else if(sph_bc%iflag_cmb .eq. iflag_fixed_field &
240+
! & .or. sph_bc%iflag_cmb .eq. iflag_evolve_field) then
241+
else
242+
call sph_out_fix_scl_val_diffuse2 &
243+
& (sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
244+
& sph_bc%kr_out, sph_bc%r_CMB, sph_bc%fdm2_fix_fld_CMB, &
245+
& CMB_Sspec%S_BC, coef_diffuse, k_ratio, dk_dr, &
246+
& scl_rj, dfs_rj)
247+
end if
248+
!
249+
end subroutine sel_CMB_sph_scalar_val_diffuse
250+
!
251+
! -----------------------------------------------------------------------
252+
!
253+
end module sel_sph_exp_scl_boundaries

‎src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_diffusion.f90‎

Lines changed: 36 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -8,18 +8,14 @@
88
!!
99
!!@verbatim
1010
!! subroutine cal_sph_nod_scalar_diffuse2(kr_in, kr_out, coef_d, &
11-
!! & is_fld, is_diffuse, nidx_rj, ar_1d_rj, g_sph_rj, &
12-
!! & d1nod_mat_fdm_2, d2nod_mat_fdm_2, &
13-
!! & n_point, ntot_phys_rj, d_rj)
11+
!! & nnod_rj, nidx_rj, ar_1d_rj, g_sph_rj, &
12+
!! & d1nod_mat_fdm_2, d2nod_mat_fdm_2, scl_rj, dfs_rj)
1413
!! subroutine cal_sph_nod_scl_val_diffuse2 &
1514
!! & (kr_in, kr_out, coef_d, k_ratio, dk_dr, &
16-
!! & is_fld, is_diffuse, nidx_rj, ar_1d_rj, g_sph_rj, &
17-
!! & d1nod_mat_fdm_2, d2nod_mat_fdm_2, &
18-
!! & n_point, ntot_phys_rj, d_rj)
15+
!! & nnod_rj, nidx_rj, ar_1d_rj, g_sph_rj, &
16+
!! & d1nod_mat_fdm_2, d2nod_mat_fdm_2, scl_rj, dfs_rj)
1917
!! integer(kind = kint), intent(in) :: kr_in, kr_out
20-
!! integer(kind = kint), intent(in) :: is_fld
21-
!! integer(kind = kint), intent(in) :: is_diffuse
22-
!! integer(kind = kint), intent(in) :: n_point, ntot_phys_rj
18+
!! integer(kind = kint), intent(in) :: nnod_rj
2319
!! integer(kind = kint), intent(in) :: nidx_rj(2)
2420
!! real(kind = kreal), intent(in) :: coef_d
2521
!! real(kind = kreal), intent(in) :: k_ratio(0:nidx_rj(1))
@@ -30,7 +26,8 @@
3026
!! & :: d1nod_mat_fdm_2(nidx_rj(1),-1:1)
3127
!! real(kind = kreal), intent(in) &
3228
!! & :: d2nod_mat_fdm_2(nidx_rj(1),-1:1)
33-
!! real(kind = kreal), intent(inout) :: d_rj(n_point,ntot_phys_rj)
29+
!! real(kind = kreal), intent(in) :: scl_rj(nnod_rj)
30+
!! real(kind = kreal), intent(inout) :: dfs_rj(nnod_rj)
3431
!!
3532
!! subroutine cal_sph_nod_vect_diffuse2(kr_in, kr_out, coef_d, &
3633
!! & is_fld, is_diffuse, nidx_rj, ar_1d_rj, g_sph_rj, &
@@ -61,14 +58,11 @@ module cal_sph_exp_diffusion
6158
! -----------------------------------------------------------------------
6259
!
6360
subroutine cal_sph_nod_scalar_diffuse2(kr_in, kr_out, coef_d, &
64-
& is_fld, is_diffuse, nidx_rj, ar_1d_rj, g_sph_rj, &
65-
& d1nod_mat_fdm_2, d2nod_mat_fdm_2, &
66-
& n_point, ntot_phys_rj, d_rj)
61+
& nnod_rj, nidx_rj, ar_1d_rj, g_sph_rj, &
62+
& d1nod_mat_fdm_2, d2nod_mat_fdm_2, scl_rj, dfs_rj)
6763
!
6864
integer(kind = kint), intent(in) :: kr_in, kr_out
69-
integer(kind = kint), intent(in) :: is_fld
70-
integer(kind = kint), intent(in) :: is_diffuse
71-
integer(kind = kint), intent(in) :: n_point, ntot_phys_rj
65+
integer(kind = kint), intent(in) :: nnod_rj
7266
integer(kind = kint), intent(in) :: nidx_rj(2)
7367
real(kind = kreal), intent(in) :: coef_d
7468
real(kind = kreal), intent(in) :: ar_1d_rj(nidx_rj(1),3)
@@ -77,8 +71,9 @@ subroutine cal_sph_nod_scalar_diffuse2(kr_in, kr_out, coef_d, &
7771
& :: d1nod_mat_fdm_2(nidx_rj(1),-1:1)
7872
real(kind = kreal), intent(in) &
7973
& :: d2nod_mat_fdm_2(nidx_rj(1),-1:1)
74+
real(kind = kreal), intent(in) :: scl_rj(nnod_rj)
8075
!
81-
real(kind = kreal), intent(inout) :: d_rj(n_point,ntot_phys_rj)
76+
real(kind = kreal), intent(inout) :: dfs_rj(nnod_rj)
8277
!
8378
real(kind = kreal) :: d1s_dr1
8479
real(kind = kreal) :: d2s_dr2
@@ -96,16 +91,15 @@ subroutine cal_sph_nod_scalar_diffuse2(kr_in, kr_out, coef_d, &
9691
j = mod((inod-1),nidx_rj(2)) + 1
9792
k = 1 + (inod- j) / nidx_rj(2)
9893
!
99-
d1s_dr1 = d1nod_mat_fdm_2(k,-1) * d_rj(i_n1,is_fld ) &
100-
& + d1nod_mat_fdm_2(k, 0) * d_rj(inod,is_fld ) &
101-
& + d1nod_mat_fdm_2(k, 1) * d_rj(i_p1,is_fld )
102-
d2s_dr2 = d2nod_mat_fdm_2(k,-1) * d_rj(i_n1,is_fld ) &
103-
& + d2nod_mat_fdm_2(k, 0) * d_rj(inod,is_fld ) &
104-
& + d2nod_mat_fdm_2(k, 1) * d_rj(i_p1,is_fld )
94+
d1s_dr1 = d1nod_mat_fdm_2(k,-1) * scl_rj(i_n1) &
95+
& + d1nod_mat_fdm_2(k, 0) * scl_rj(inod) &
96+
& + d1nod_mat_fdm_2(k, 1) * scl_rj(i_p1)
97+
d2s_dr2 = d2nod_mat_fdm_2(k,-1) * scl_rj(i_n1) &
98+
& + d2nod_mat_fdm_2(k, 0) * scl_rj(inod) &
99+
& + d2nod_mat_fdm_2(k, 1) * scl_rj(i_p1)
105100
!
106-
d_rj(inod,is_diffuse ) &
107-
& = coef_d * (d2s_dr2 + two*ar_1d_rj(k,1) * d1s_dr1 &
108-
& - g_sph_rj(j,3)*ar_1d_rj(k,2)*d_rj(inod,is_fld) )
101+
dfs_rj(inod) = coef_d * (d2s_dr2 + two*ar_1d_rj(k,1) * d1s_dr1 &
102+
& - g_sph_rj(j,3)*ar_1d_rj(k,2)*scl_rj(inod))
109103
end do
110104
!$omp end parallel do
111105
!
@@ -115,14 +109,11 @@ end subroutine cal_sph_nod_scalar_diffuse2
115109
!
116110
subroutine cal_sph_nod_scl_val_diffuse2 &
117111
& (kr_in, kr_out, coef_d, k_ratio, dk_dr, &
118-
& is_fld, is_diffuse, nidx_rj, ar_1d_rj, g_sph_rj, &
119-
& d1nod_mat_fdm_2, d2nod_mat_fdm_2, &
120-
& n_point, ntot_phys_rj, d_rj)
112+
& nnod_rj, nidx_rj, ar_1d_rj, g_sph_rj, &
113+
& d1nod_mat_fdm_2, d2nod_mat_fdm_2, scl_rj, dfs_rj)
121114
!
122115
integer(kind = kint), intent(in) :: kr_in, kr_out
123-
integer(kind = kint), intent(in) :: is_fld
124-
integer(kind = kint), intent(in) :: is_diffuse
125-
integer(kind = kint), intent(in) :: n_point, ntot_phys_rj
116+
integer(kind = kint), intent(in) :: nnod_rj
126117
integer(kind = kint), intent(in) :: nidx_rj(2)
127118
!
128119
real(kind = kreal), intent(in) :: coef_d
@@ -134,8 +125,9 @@ subroutine cal_sph_nod_scl_val_diffuse2 &
134125
& :: d1nod_mat_fdm_2(nidx_rj(1),-1:1)
135126
real(kind = kreal), intent(in) &
136127
& :: d2nod_mat_fdm_2(nidx_rj(1),-1:1)
128+
real(kind = kreal), intent(in) :: scl_rj(nnod_rj)
137129
!
138-
real(kind = kreal), intent(inout) :: d_rj(n_point,ntot_phys_rj)
130+
real(kind = kreal), intent(inout) :: dfs_rj(nnod_rj)
139131
!
140132
real(kind = kreal) :: d1s_dr1
141133
real(kind = kreal) :: d2s_dr2
@@ -153,17 +145,17 @@ subroutine cal_sph_nod_scl_val_diffuse2 &
153145
j = mod((inod-1),nidx_rj(2)) + 1
154146
k = 1 + (inod- j) / nidx_rj(2)
155147
!
156-
d1s_dr1 = d1nod_mat_fdm_2(k,-1) * d_rj(i_n1,is_fld ) &
157-
& + d1nod_mat_fdm_2(k, 0) * d_rj(inod,is_fld ) &
158-
& + d1nod_mat_fdm_2(k, 1) * d_rj(i_p1,is_fld )
159-
d2s_dr2 = d2nod_mat_fdm_2(k,-1) * d_rj(i_n1,is_fld ) &
160-
& + d2nod_mat_fdm_2(k, 0) * d_rj(inod,is_fld ) &
161-
& + d2nod_mat_fdm_2(k, 1) * d_rj(i_p1,is_fld )
162-
!
163-
d_rj(inod,is_diffuse) &
164-
& = coef_d*k_ratio(k) * (d2s_dr2 + two*ar_1d_rj(k,1) * d1s_dr1 &
165-
& - g_sph_rj(j,3)*ar_1d_rj(k,2)*d_rj(inod,is_fld)) &
166-
& + coef_d * dk_dr(k) * d1s_dr1
148+
d1s_dr1 = d1nod_mat_fdm_2(k,-1) * scl_rj(i_n1) &
149+
& + d1nod_mat_fdm_2(k, 0) * scl_rj(inod) &
150+
& + d1nod_mat_fdm_2(k, 1) * scl_rj(i_p1)
151+
d2s_dr2 = d2nod_mat_fdm_2(k,-1) * scl_rj(i_n1) &
152+
& + d2nod_mat_fdm_2(k, 0) * scl_rj(inod) &
153+
& + d2nod_mat_fdm_2(k, 1) * scl_rj(i_p1)
154+
!
155+
dfs_rj(inod) = coef_d*k_ratio(k) &
156+
& * (d2s_dr2 + two*ar_1d_rj(k,1) * d1s_dr1 &
157+
& - g_sph_rj(j,3) * ar_1d_rj(k,2) * scl_rj(inod)) &
158+
& + coef_d * dk_dr(k) * d1s_dr1
167159
end do
168160
!$omp end parallel do
169161
!

0 commit comments

Comments
 (0)