@@ -79,22 +79,22 @@ contains
7979 real (wp), optional, intent (in ) :: div
8080
8181 integer :: i, j
82- integer :: m, n
82+ integer :: local_m, local_n
8383 real (wp) :: c
8484
85- m = size (A, 1 )
86- n = size (A, 2 )
85+ local_m = size (A, 1 )
86+ local_n = size (A, 2 )
8787
8888 if (present (div)) then
8989 c = div
9090 else
9191 c = 1._wp
9292 end if
9393
94- print * , m, n
94+ print * , local_m, local_n
9595
96- do i = 1 , m
97- do j = 1 , n
96+ do i = 1 , local_m
97+ do j = 1 , local_n
9898 write (* , fmt= " (F12.4)" , advance= " no" ) A(i, j)/ c
9999 end do
100100 write (* , fmt= " (A1)" ) " "
@@ -251,10 +251,10 @@ contains
251251 end subroutine s_int_to_str
252252
253253 !> Computes the Simpson weights for quadrature
254- subroutine s_simpson (weight , R0 )
254+ subroutine s_simpson (local_weight , local_R0 )
255255
256- real (wp), dimension (:), intent (inout ) :: weight
257- real (wp), dimension (:), intent (inout ) :: R0
256+ real (wp), dimension (:), intent (inout ) :: local_weight
257+ real (wp), dimension (:), intent (inout ) :: local_R0
258258
259259 integer :: ir
260260 real (wp) :: R0mn, R0mx, dphi, tmp, sd
@@ -268,7 +268,7 @@ contains
268268 do ir = 1 , nb
269269 phi(ir) = log (R0mn) &
270270 + (ir - 1._wp )* log (R0mx/ R0mn)/ (nb - 1._wp )
271- R0 (ir) = exp (phi(ir))
271+ local_R0 (ir) = exp (phi(ir))
272272 end do
273273 dphi = phi(2 ) - phi(1 )
274274
@@ -277,15 +277,15 @@ contains
277277 ! Gaussian
278278 tmp = exp(-0.5_wp*(phi(ir)/sd)**2)/sqrt(2._wp*pi)/sd
279279 if (mod(ir, 2) == 0) then
280- weight (ir) = tmp*4._wp*dphi/3._wp
280+ local_weight (ir) = tmp*4._wp*dphi/3._wp
281281 else
282- weight (ir) = tmp*2._wp*dphi/3._wp
282+ local_weight (ir) = tmp*2._wp*dphi/3._wp
283283 end if
284284 end do
285285 tmp = exp(-0.5_wp*(phi(1)/sd)**2)/sqrt(2._wp*pi)/sd
286- weight (1) = tmp*dphi/3._wp
286+ local_weight (1) = tmp*dphi/3._wp
287287 tmp = exp(-0.5_wp*(phi(nb)/sd)**2)/sqrt(2._wp*pi)/sd
288- weight (nb) = tmp*dphi/3._wp
288+ local_weight (nb) = tmp*dphi/3._wp
289289 end subroutine s_simpson
290290
291291 !> This procedure computes the cross product of two vectors.
@@ -318,40 +318,40 @@ contains
318318 !> This procedure creates a transformation matrix.
319319 !! @param p Parameters for the transformation.
320320 !! @return Transformation matrix.
321- pure function f_create_transform_matrix(p , center) result(out_matrix)
321+ pure function f_create_transform_matrix(param , center) result(out_matrix)
322322
323- type(ic_model_parameters), intent(in) :: p
323+ type(ic_model_parameters), intent(in) :: param
324324 real(wp), dimension(1:3), optional, intent(in) :: center
325325 real(wp), dimension(1:4, 1:4) :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix
326326
327327 sc = transpose(reshape([ &
328- p %scale(1), 0._wp, 0._wp, 0._wp, &
329- 0._wp, p %scale(2), 0._wp, 0._wp, &
330- 0._wp, 0._wp, p %scale(3), 0._wp, &
328+ param %scale(1), 0._wp, 0._wp, 0._wp, &
329+ 0._wp, param %scale(2), 0._wp, 0._wp, &
330+ 0._wp, 0._wp, param %scale(3), 0._wp, &
331331 0._wp, 0._wp, 0._wp, 1._wp], shape(sc)))
332332
333333 rz = transpose(reshape([ &
334- cos(p %rotate(3)), -sin(p %rotate(3)), 0._wp, 0._wp, &
335- sin(p %rotate(3)), cos(p %rotate(3)), 0._wp, 0._wp, &
334+ cos(param %rotate(3)), -sin(param %rotate(3)), 0._wp, 0._wp, &
335+ sin(param %rotate(3)), cos(param %rotate(3)), 0._wp, 0._wp, &
336336 0._wp, 0._wp, 1._wp, 0._wp, &
337337 0._wp, 0._wp, 0._wp, 1._wp], shape(rz)))
338338
339339 rx = transpose(reshape([ &
340340 1._wp, 0._wp, 0._wp, 0._wp, &
341- 0._wp, cos(p %rotate(1)), -sin(p %rotate(1)), 0._wp, &
342- 0._wp, sin(p %rotate(1)), cos(p %rotate(1)), 0._wp, &
341+ 0._wp, cos(param %rotate(1)), -sin(param %rotate(1)), 0._wp, &
342+ 0._wp, sin(param %rotate(1)), cos(param %rotate(1)), 0._wp, &
343343 0._wp, 0._wp, 0._wp, 1._wp], shape(rx)))
344344
345345 ry = transpose(reshape([ &
346- cos(p %rotate(2)), 0._wp, sin(p %rotate(2)), 0._wp, &
346+ cos(param %rotate(2)), 0._wp, sin(param %rotate(2)), 0._wp, &
347347 0._wp, 1._wp, 0._wp, 0._wp, &
348- -sin(p %rotate(2)), 0._wp, cos(p %rotate(2)), 0._wp, &
348+ -sin(param %rotate(2)), 0._wp, cos(param %rotate(2)), 0._wp, &
349349 0._wp, 0._wp, 0._wp, 1._wp], shape(ry)))
350350
351351 tr = transpose(reshape([ &
352- 1._wp, 0._wp, 0._wp, p %translate(1), &
353- 0._wp, 1._wp, 0._wp, p %translate(2), &
354- 0._wp, 0._wp, 1._wp, p %translate(3), &
352+ 1._wp, 0._wp, 0._wp, param %translate(1), &
353+ 0._wp, 1._wp, 0._wp, param %translate(2), &
354+ 0._wp, 0._wp, 1._wp, param %translate(3), &
355355 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
356356
357357 if (present(center)) then
@@ -484,18 +484,18 @@ contains
484484 !! @param x is the input value
485485 !! @param l is the degree
486486 !! @return P is the unassociated legendre polynomial evaluated at x
487- pure recursive function unassociated_legendre(x, l) result(P )
487+ pure recursive function unassociated_legendre(x, l) result(result_P )
488488
489489 integer, intent(in) :: l
490490 real(wp), intent(in) :: x
491- real(wp) :: P
491+ real(wp) :: result_P
492492
493493 if (l == 0) then
494- P = 1._wp
494+ result_P = 1._wp
495495 else if (l == 1) then
496- P = x
496+ result_P = x
497497 else
498- P = ((2*l - 1)*x*unassociated_legendre(x, l - 1) - (l - 1)*unassociated_legendre(x, l - 2))/l
498+ result_P = ((2*l - 1)*x*unassociated_legendre(x, l - 1) - (l - 1)*unassociated_legendre(x, l - 2))/l
499499 end if
500500
501501 end function unassociated_legendre
@@ -504,20 +504,20 @@ contains
504504 !! @param x is the x coordinate
505505 !! @param phi is the phi coordinate
506506 !! @param l is the degree
507- !! @param m is the order
507+ !! @param m_order is the order
508508 !! @return Y is the spherical harmonic function evaluated at x and phi
509- pure recursive function spherical_harmonic_func(x, phi, l, m ) result(Y)
509+ pure recursive function spherical_harmonic_func(x, phi, l, m_order ) result(Y)
510510
511- integer, intent(in) :: l, m
511+ integer, intent(in) :: l, m_order
512512 real(wp), intent(in) :: x, phi
513- real(wp) :: Y, prefactor, pi
514-
515- pi = acos(-1._wp)
516- prefactor = sqrt((2*l + 1)/(4*pi )*factorial(l - m )/factorial(l + m ));
517- if (m == 0) then
518- Y = prefactor*associated_legendre(x, l, m );
519- elseif (m > 0) then
520- Y = (-1._wp)**m *sqrt(2._wp)*prefactor*associated_legendre(x, l, m )*cos(m *phi);
513+ real(wp) :: Y, prefactor, local_pi
514+
515+ local_pi = acos(-1._wp)
516+ prefactor = sqrt((2*l + 1)/(4*local_pi )*factorial(l - m_order )/factorial(l + m_order ));
517+ if (m_order == 0) then
518+ Y = prefactor*associated_legendre(x, l, m_order );
519+ elseif (m_order > 0) then
520+ Y = (-1._wp)**m_order *sqrt(2._wp)*prefactor*associated_legendre(x, l, m_order )*cos(m_order *phi);
521521 end if
522522
523523 end function spherical_harmonic_func
@@ -528,54 +528,54 @@ contains
528528 !! @param l is the degree
529529 !! @param m is the order
530530 !! @return P is the associated legendre polynomial evaluated at x
531- pure recursive function associated_legendre(x, l, m ) result(P )
531+ pure recursive function associated_legendre(x, l, m_order ) result(result_P )
532532
533- integer, intent(in) :: l, m
533+ integer, intent(in) :: l, m_order
534534 real(wp), intent(in) :: x
535- real(wp) :: P
536-
537- if (m <= 0 .and. l <= 0) then
538- P = 1;
539- elseif (l == 1 .and. m <= 0) then
540- P = x;
541- elseif (l == 1 .and. m == 1) then
542- P = -(1 - x**2)**(1._wp/2._wp);
543- elseif (m == l) then
544- P = (-1)**l*double_factorial(2*l - 1)*(1 - x**2)**(l/2);
545- elseif (m == l - 1) then
546- P = x*(2*l - 1)*associated_legendre(x, l - 1, l - 1);
535+ real(wp) :: result_P
536+
537+ if (m_order <= 0 .and. l <= 0) then
538+ result_P = 1;
539+ elseif (l == 1 .and. m_order <= 0) then
540+ result_P = x;
541+ elseif (l == 1 .and. m_order == 1) then
542+ result_P = -(1 - x**2)**(1._wp/2._wp);
543+ elseif (m_order == l) then
544+ result_P = (-1)**l*double_factorial(2*l - 1)*(1 - x**2)**(l/2);
545+ elseif (m_order == l - 1) then
546+ result_P = x*(2*l - 1)*associated_legendre(x, l - 1, l - 1);
547547 else
548- P = ((2*l - 1)*x*associated_legendre(x, l - 1, m ) - (l + m - 1)*associated_legendre(x, l - 2, m ))/(l - m );
548+ result_P = ((2*l - 1)*x*associated_legendre(x, l - 1, m_order ) - (l + m_order - 1)*associated_legendre(x, l - 2, m_order ))/(l - m_order );
549549 end if
550550
551551 end function associated_legendre
552552
553553 !> This function calculates the double factorial value of an integer
554- !! @param n is the input integer
554+ !! @param n_in is the input integer
555555 !! @return R is the double factorial value of n
556- pure elemental function double_factorial(n ) result(R )
556+ pure elemental function double_factorial(n_in ) result(R_result )
557557
558- integer, intent(in) :: n
558+ integer, intent(in) :: n_in
559559 integer, parameter :: int64_kind = selected_int_kind(18) ! 18 bytes for 64-bit integer
560- integer(kind=int64_kind) :: R
560+ integer(kind=int64_kind) :: R_result
561561 integer :: i
562562
563- R = product((/(i, i=n , 1, -2)/))
563+ R_result = product((/(i, i=n_in , 1, -2)/))
564564
565565 end function double_factorial
566566
567567 !> The following function calculates the factorial value of an integer
568- !! @param n is the input integer
568+ !! @param n_in is the input integer
569569 !! @return R is the factorial value of n
570- pure elemental function factorial(n ) result(R )
570+ pure elemental function factorial(n_in ) result(R_result )
571571
572- integer, intent(in) :: n
572+ integer, intent(in) :: n_in
573573 integer, parameter :: int64_kind = selected_int_kind(18) ! 18 bytes for 64-bit integer
574- integer(kind=int64_kind) :: R
574+ integer(kind=int64_kind) :: R_result
575575
576576 integer :: i
577577
578- R = product((/(i, i=n , 1, -1)/))
578+ R_result = product((/(i, i=n_in , 1, -1)/))
579579
580580 end function factorial
581581
0 commit comments