53 real(4) :: lf = 0.0, uf = 0.0
56 real(8) :: llc_cart = 0.0
58 real(8) :: llc_nat = 0.0
61 real(8),
dimension(:),
pointer :: r => null()
62 real(8),
dimension(:),
pointer :: h => null()
64 real,
dimension(:),
pointer :: ddn => null()
65 real,
dimension(:),
pointer :: dup => null()
66 real,
dimension(:),
pointer :: dsup => null()
67 real,
dimension(:),
pointer :: dsdn => null()
68 real,
dimension(:),
pointer :: b_m => null()
69 real,
dimension(:),
pointer :: b_mdn => null()
70 real,
dimension(:),
pointer :: b_1d => null()
71 real,
dimension(:),
pointer :: b_didup => null()
72 real,
dimension(:),
pointer :: b_diddn => null()
73 real,
dimension(:,:),
pointer:: zlincoeff => null()
74 real,
dimension(:,:),
pointer:: dzlincoeff => null()
75 real,
dimension(:,:),
pointer:: zuincoeff => null()
76 real,
dimension(:,:),
pointer:: dzuincoeff => null()
78 real :: b_ad, b_bd, b_cd
83 integer(kind=4) :: b_r, b_rb
85 character(len=16)::
type =
'' 86 real(8),
allocatable :: zeta
87 real(8),
allocatable :: xi
92 logical:: lower_boundary = .false.
93 logical:: upper_boundary = .false.
94 logical:: no_mans_land = .false.
96 procedure(),
pointer,
nopass :: currenttocartesian => null()
98 procedure(func_interface),
pointer :: h2c => null()
99 procedure(func_interface),
pointer :: h2f => null()
100 procedure(func_interface),
pointer :: h31c => null()
101 procedure(func_interface),
pointer :: h31f => null()
102 procedure(func_interface),
pointer :: h32f => null()
103 procedure(func_interface),
pointer :: h32c => null()
104 procedure(func_interface),
pointer :: dvol1c => null()
105 procedure(func_interface),
pointer :: dvol1f => null()
106 procedure(func_interface),
pointer :: dvol2c => null()
107 procedure(func_interface),
pointer :: dvol2f => null()
108 procedure(func_interface),
pointer :: dvol3c => null()
109 procedure(func_interface),
pointer :: dvol3f => null()
110 procedure(func_interface),
pointer :: dar1c => null()
111 procedure(func_interface),
pointer :: dar1f => null()
112 procedure(func_interface),
pointer :: dar2c => null()
113 procedure(func_interface),
pointer :: dar2f => null()
114 procedure(func_interface),
pointer :: dar31c => null()
115 procedure(func_interface),
pointer :: dar31f => null()
116 procedure(func_interface),
pointer :: dar32c => null()
117 procedure(func_interface),
pointer :: dar32f => null()
119 procedure :: init => init_mesh
120 procedure :: deallocate => deallocate_coord
121 procedure :: print => print_mesh
122 procedure(init_geometryinterface),
deferred :: init_geometry
125 subroutine init_geometryinterface(self, idir, origin)
128 integer,
intent(in) :: idir
129 real(8),
intent(in) :: origin
132 function func_interface(self, i)
result(x)
135 integer,
intent(in) :: i
144 integer :: cartesian = 1
145 integer :: spherical = 2
146 integer :: cylindrical = 3
148 type(mesh_type),
public :: mesh_types
155 procedure :: init_geometry => init_cartesian_geometry
166 procedure :: init_geometry => init_cylindrical_geometry
180 procedure :: init_geometry => init_spherical_geometry
184 module procedure :: cylindricaltocartesian4, cylindricaltocartesian8
186 interface cylindricaltospherical
187 module procedure :: cylindricaltospherical4, cylindricaltospherical8
189 interface sphericaltocartesian
190 module procedure :: sphericaltocartesian4, sphericaltocartesian8
192 interface cartesiantospherical
193 module procedure :: cartesiantospherical4, cartesiantospherical8
195 interface cartesiantocylindrical
196 module procedure :: cartesiantocylindrical4, cartesiantocylindrical8
198 interface sphericaltocylindrical
199 module procedure :: sphericaltocylindrical4, sphericaltocylindrical8
201 interface currenttocartesian
202 module procedure :: currenttocartesian4, currenttocartesian8
204 interface cartesiantocurrent
205 module procedure :: cartesiantocurrent4, cartesiantocurrent8
208 module procedure :: round4, round8
211 real(8),
parameter :: pi = 3.14159265358979323846
213 PUBLIC meshfactory, meshrecycler, &
215 cylindricaltospherical, sphericaltocartesian, sphericaltocylindrical, &
216 currenttocartesian, cartesiantocurrent, round
224 SUBROUTINE meshfactory (meshtype, collection, n, ng, s, p, llc_native, nml)
225 class(
mesh_t),
dimension(:),
pointer :: collection
226 integer,
intent(in) :: meshtype
227 integer,
dimension(3),
optional :: ng
228 integer,
dimension(3) :: ngi
229 integer,
dimension(3),
optional :: n
230 integer,
dimension(3) :: ni
231 real(8),
dimension(3),
optional :: s, p
232 real(8),
dimension(3) :: si, pi
233 logical,
optional :: nml
238 real(8),
dimension(3),
optional :: llc_native
239 real(8),
dimension(3) :: llc_native_i
257 if (
present(llc_native))
then 258 llc_native_i = llc_native
260 llc_native_i(:) = 0.0
262 if (
present(ng))
then 269 if (
associated(collection))
call meshrecycler(collection)
273 if (meshtype .eq. mesh_types%Cartesian)
then 275 collection%type =
'Cartesian_mesh' 276 else if ( meshtype .eq. mesh_types%spherical )
then 278 collection%type =
'spherical_mesh' 279 else if ( meshtype .eq. mesh_types%cylindrical )
then 281 collection%type =
'cylindrical_mesh' 283 call mpi%abort(
"Invalid mesh type. Abort!")
285 if (
present(nml))
then 286 collection%no_mans_land = nml
291 collection(i)%idir = i
292 call collection(i)%init(ni(i),ngi(i),si(i),pi(i),llc_native_i(i))
295 call collection(i)%init_geometry(i,llc_native_i(i))
298 END SUBROUTINE meshfactory
303 SUBROUTINE meshrecycler (collection)
304 class(
mesh_t),
dimension(:),
pointer :: collection
307 do i=1,
size(collection)
308 call collection(i)%deallocate
310 deallocate(collection)
312 END SUBROUTINE meshrecycler
319 SUBROUTINE init_mesh (self, n, ng, s, p, llc_native)
322 integer,
optional :: n, ng
323 real(8),
optional :: s, p
324 real(8),
optional :: llc_native
327 call trace_begin (
'mesh%init')
332 if (
present(ng))
then 341 if (
present(llc_native))
then 342 self%llc_nat = llc_native
347 self%li = self%lb + self%ng
348 self%lo = self%li - 1
349 self%ui = self%lo + self%n
350 self%uo = self%ui + 1
351 self%ub = self%ui + self%ng
353 if (self%n == 1)
then 359 allocate (self%r(self%gn))
361 self%o = (self%li+self%ui)*0.5
364 self%centre_nat = volumecentroid(self)
369 if (self%no_mans_land)
then 371 if (self%n == 1) h = 1.0_8
372 self%nc = self%ui - self%li + 1
375 self%nc = max(1, self%ui - self%li)
377 self%d = self%s/max(1,self%nc)
379 self%r(i) =
real(i - h - self%ng,kind=8) * self%d + self%llc_nat - self%centre_nat
381 self%lf = -self%s*0.5_8
382 self%uf = +self%s*0.5_8
384 END SUBROUTINE init_mesh
389 SUBROUTINE deallocate_coord (self)
392 if (
associated(self%r))
deallocate (self%r)
393 if (
associated(self%h))
deallocate (self%h)
394 if (
associated(self%ddn))
deallocate (self%ddn)
395 if (
associated(self%dup))
deallocate (self%dup)
396 END SUBROUTINE deallocate_coord
401 SUBROUTINE print_mesh (self)
406 if (io%verbose <= 0)
return 407 print*,
'::mesh%print::' 410 print*,
'Cartesian', self%idir
412 print*,
'spherical', self%idir
414 print*,
'cylindrical', self%idir
416 print
'(a,3i8)',
"n, ng, gn: ", self%n, self%ng, self%gn
417 print
'(a,6i8)',
"lb, lo, li, ub, uo, ui: ", self%lb, self%lo, self%li, self%ub, self%uo, self%ui
418 print
'(a,3x,3g16.8)',
"lf, uf, p: ",self%lf, self%uf, self%p
419 print
'(a,3x,3g16.8)',
"s, d, b: ", self%s, self%d, self%b
420 print
'(a,3x,4g16.8)',
"llc_nat, centre_nat, llc_cart, o: ",self%llc_nat, self%centre_nat, self%llc_cart, self%o
422 print
'(i4,f14.8)',i,self%r(i)
425 if (
allocated(self%zeta)) print*,
"phi axis offset: ", self%zeta
426 if (
allocated(self%xi)) print*,
"theta axis offset: ", self%xi
429 END SUBROUTINE print_mesh
435 SUBROUTINE init_cartesian_geometry (self, idir, origin)
437 integer,
intent(in) :: idir
438 real(8),
intent(in) :: origin
448 self%h31c => func_one
449 self%h31f => func_one
450 self%dvol1c => func_dx
451 self%dvol1f => func_dx
453 self%dar1c => func_one
454 self%dar1f => func_one
455 self%dar2c => func_one
456 self%dar2f => func_one
457 self%dar31c => func_one
458 self%dar31f => func_one
460 self%h32c => func_one
461 self%h32f => func_one
462 self%dvol2c => func_dx
463 self%dvol2f => func_dx
465 self%dar32c => func_one
466 self%dar32f => func_one
468 self%dvol3c => func_dx
469 self%dvol3f => func_dx
472 END SUBROUTINE init_cartesian_geometry
478 SUBROUTINE init_cylindrical_geometry (self, idir, origin)
480 integer,
intent(in) :: idir
481 real(8),
intent(in) :: origin
487 self%h31c => func_one
488 self%h31f => func_one
489 self%dvol1c => func_dx
490 self%dvol1f => func_dx
492 self%dar1c => func_one
493 self%dar1f => func_one
494 self%dar2c => func_one
495 self%dar2f => func_one
496 self%dar31c => func_one
497 self%dar31f => func_one
499 self%h32c => func_radius_c
500 self%h32f => func_radius_f
501 self%dvol2c => func_dr2_c
502 self%dvol2f => func_dr2_f
504 self%dar32c => func_dar32_c
505 self%dar32f => func_dar32_f
507 self%dvol3c => func_dx
508 self%dvol3f => func_dx
515 END SUBROUTINE init_cylindrical_geometry
521 SUBROUTINE init_spherical_geometry (self, idir, origin)
523 integer,
intent(in) :: idir
524 real(8),
intent(in) :: origin
529 self%h2c => func_radius_c
530 self%h2f => func_radius_f
531 self%h31c => func_radius_c
532 self%h31f => func_radius_f
533 self%dvol1c => func_dr3_c
534 self%dvol1f => func_dr3_f
536 self%dar1c => func_dar1_c
537 self%dar1f => func_dar1_f
538 self%dar2c => func_dar2_c
539 self%dar2f => func_dar2_f
540 self%dar31c => func_dar31_c
541 self%dar31f => func_dar31_f
543 self%h32c => func_sintheta_c
544 self%h32f => func_sintheta_f
545 self%dvol2c => func_dcostheta_c
546 self%dvol2f => func_dcostheta_f
548 self%dar32c => func_dar32_c
549 self%dar32f => func_dar32_f
555 self%dvol3c => func_dx
556 self%dvol3f => func_dx
563 END SUBROUTINE init_spherical_geometry
568 FUNCTION func_one (self, i)
RESULT(one)
570 integer,
intent(in) :: i
575 END FUNCTION func_one
577 FUNCTION func_dx (self, i)
RESULT(dx)
579 integer,
intent(in) :: i
586 FUNCTION func_position_c (self, i)
RESULT(r)
588 integer,
intent(in) :: i
591 r = self%centre_nat + self%r(i)
593 END FUNCTION func_position_c
595 FUNCTION func_position_f (self, i)
RESULT(r)
597 integer,
intent(in) :: i
600 if (i /= self%gn .or. self%gn == 1)
then 601 r = self%centre_nat + self%r(i) - 0.5 * self%d
603 r = self%centre_nat + self%r(i-1) + 0.5 * self%d
606 END FUNCTION func_position_f
608 FUNCTION func_radius_c (self, i)
RESULT(r)
610 integer,
intent(in) :: i
613 r = abs(func_position_c(self, i))
615 END FUNCTION func_radius_c
617 FUNCTION func_radius_f (self, i)
RESULT(r)
619 integer,
intent(in) :: i
622 r = abs(func_position_f(self, i))
624 END FUNCTION func_radius_f
626 FUNCTION func_sintheta_c (self, j)
RESULT(sintc)
628 integer,
intent(in) :: j
631 sintc = abs(sin(func_position_c(self,j)))
633 END FUNCTION func_sintheta_c
635 FUNCTION func_sintheta_f (self, j)
RESULT(sintf)
637 integer,
intent(in) :: j
640 sintf = abs(sin(func_position_f(self,j)))
642 END FUNCTION func_sintheta_f
651 FUNCTION func_dr2_c (self, j)
RESULT(dvol2c)
653 integer,
intent(in) :: j
656 dvol2c = 0.5 * abs( self%h32f(j+1) * func_position_f(self,j+1) &
657 - self%h32f(j ) * func_position_f(self,j ))
659 END FUNCTION func_dr2_c
661 FUNCTION func_dr2_f (self, j)
RESULT(dvol2f)
663 integer,
intent(in) :: j
666 dvol2f = 0.5 * abs( self%h32c(j ) * func_position_c(self,j ) &
667 - self%h32c(j-1) * func_position_c(self,j-1))
669 END FUNCTION func_dr2_f
671 FUNCTION func_dr3_c (self, i)
RESULT(dvol1c)
673 integer,
intent(in) :: i
676 dvol1c = self%h2f(i+1) * self%h31f(i+1) * func_position_f(self,i+1) / 3.0d0 &
677 - self%h2f(i ) * self%h31f(i ) * func_position_f(self,i ) / 3.0d0
679 END FUNCTION func_dr3_c
681 FUNCTION func_dr3_f (self, i)
RESULT(dvol1f)
683 integer,
intent(in) :: i
686 dvol1f = self%h2c(i ) * self%h31c(i ) * func_position_c(self,i ) / 3.0d0 &
687 - self%h2c(i-1) * self%h31c(i-1) * func_position_c(self,i-1) / 3.0d0
689 END FUNCTION func_dr3_f
691 FUNCTION func_dcostheta_c (self, j)
RESULT(dvol2c)
693 integer,
intent(in) :: j
696 dvol2c = -cos(func_position_f(self,j+1)) + cos(func_position_f(self,j))
698 END FUNCTION func_dcostheta_c
700 FUNCTION func_dcostheta_f (self, j)
RESULT(dvol2f)
702 integer,
intent(in) :: j
705 dvol2f = -cos(func_position_c(self,j)) + cos(func_position_c(self,j-1))
706 END FUNCTION func_dcostheta_f
711 FUNCTION func_dar1_c (self, i)
RESULT(dar1c)
713 integer,
intent(in) :: i
716 dar1c = self%h2c(i) * self%h31c(i)
718 END FUNCTION func_dar1_c
720 FUNCTION func_dar1_f (self, i)
RESULT(dar1f)
722 integer,
intent(in) :: i
725 dar1f = self%h2f(i) * self%h31f(i)
727 END FUNCTION func_dar1_f
729 FUNCTION func_dar2_c (self,i)
RESULT(dar2c)
731 integer,
intent(in) :: i
734 dar2c = self%h31c(i) * self%d / self%dvol1c(i)
736 END FUNCTION func_dar2_c
738 FUNCTION func_dar2_f (self,i)
RESULT(dar2f)
740 integer,
intent(in) :: i
743 dar2f = self%h31f(i) * self%d / self%dvol1f(i)
745 END FUNCTION func_dar2_f
747 FUNCTION func_dar31_c (self,i)
RESULT(dar31c)
749 integer,
intent(in) :: i
752 dar31c = self%h2c(i) * self%d / self%dvol1c(i)
754 END FUNCTION func_dar31_c
756 FUNCTION func_dar31_f (self,i)
RESULT(dar31f)
758 integer,
intent(in) :: i
761 dar31f = self%h2f(i) * self%d / self%dvol1f(i)
763 END FUNCTION func_dar31_f
765 FUNCTION func_dar32_c (self, j)
RESULT(dar32c)
767 integer,
intent(in) :: j
770 dar32c = self%d / self%dvol2c(j)
772 END FUNCTION func_dar32_c
774 FUNCTION func_dar32_f (self, j)
RESULT(dar32f)
776 integer,
intent(in) :: j
779 dar32f = self%d / self%dvol2f(j)
781 END FUNCTION func_dar32_f
787 FUNCTION cartesiantocylindrical4 (x, y, z, OCart)
result (vec)
788 real(4),
intent(in) :: x, y, z
789 real(4),
dimension(3),
optional:: OCart
790 real(4),
dimension(3) :: vec, lOCart
793 if (
present(ocart)) locart = ocart
795 vec(1) = z - locart(3)
796 vec(2) = sqrt((x-locart(1))**2 + (y-locart(2))**2)
797 if (locart(1) == 0.0 .and. locart(2) == 0.0)
then 800 vec(3) = atan2(y-locart(2), x-locart(1))
802 if (vec(3) .lt. 0.0) vec(3) = vec(3) + 2.0 *
real(pi,kind=4)
804 END FUNCTION cartesiantocylindrical4
806 FUNCTION cartesiantocylindrical8 (x, y, z, OCart)
result (vec)
807 real(8),
intent(in) :: x, y, z
808 real(8),
dimension(3),
optional:: OCart
809 real(8),
dimension(3) :: vec, lOCart
812 if (
present(ocart)) locart = ocart
814 vec(1) = z - locart(3)
815 vec(2) = sqrt((x-locart(1))**2 + (y-locart(2))**2)
816 if (locart(1) == 0.0 .and. locart(2) == 0.0)
then 819 vec(3) = atan2(y-locart(2), x-locart(1))
821 if (vec(3) .lt. 0.0d0) vec(3) = vec(3) + 2.0d0 * pi
823 END FUNCTION cartesiantocylindrical8
825 FUNCTION cartesiantospherical4 (x, y, z, OCart)
result (vec)
826 real(4),
intent(in) :: x, y, z
827 real(4),
dimension(3),
optional:: OCart
828 real(4),
dimension(3) :: vec, lOCart
831 if (
present(ocart)) locart = ocart
833 vec(1) = sqrt((x-locart(1))**2 + (y-locart(2))**2 + (z-locart(3))**2)
834 vec(2) = pi - acos((z-locart(3)) / vec(1))
835 if (locart(1) == 0.0 .and. locart(2) == 0.0)
then 836 vec(3) = -atan2(y, x)
838 vec(3) = -atan2((y-locart(2)), (x-locart(1)))
840 if (vec(3) .lt. 0.0) vec(3) = vec(3) + 2.0 *
real(pi,kind=4)
842 END FUNCTION cartesiantospherical4
844 FUNCTION cartesiantospherical8 (x, y, z, OCart)
result (vec)
845 real(8),
intent(in) :: x, y, z
846 real(8),
dimension(3),
optional:: OCart
847 real(8),
dimension(3) :: vec, lOCart
850 if (
present(ocart)) locart = ocart
852 vec(1) = sqrt((x-locart(1))**2 + (y-locart(2))**2 + (z-locart(3))**2)
853 vec(2) = pi - acos((z-locart(3)) / vec(1))
854 if (locart(1) == 0.0 .and. locart(2) == 0.0)
then 855 vec(3) = -atan2(y, x)
857 vec(3) = -atan2((y-locart(2)), (x-locart(1)))
859 if (vec(3) .lt. 0.0d0) vec(3) = vec(3) + 2.0d0 * pi
861 END FUNCTION cartesiantospherical8
863 FUNCTION cylindricaltocartesian4 (Z, R, Phi, OCart)
result (vec)
864 real(4),
intent(in) :: Z, R, Phi
865 real(4),
dimension(3),
optional :: OCart
866 real(4),
dimension(3) :: vec, lOCart
869 if (
present(ocart)) locart = ocart
871 vec(1) = r * cos(phi) - locart(1)
872 vec(2) = r * sin(phi) - locart(2)
873 vec(3) = z - locart(3)
875 END FUNCTION cylindricaltocartesian4
877 FUNCTION cylindricaltocartesian8 (Z, R, Phi, OCart)
result (vec)
878 real(8),
intent(in) :: Z, R, Phi
879 real(8),
dimension(3),
optional :: OCart
880 real(8),
dimension(3) :: vec, lOCart
883 if (
present(ocart)) locart = ocart
885 vec(1) = r * cos(phi) - locart(1)
886 vec(2) = r * sin(phi) - locart(2)
887 vec(3) = z - locart(3)
889 END FUNCTION cylindricaltocartesian8
891 FUNCTION sphericaltocartesian4 (r, theta, phi, OCart)
result (vec)
892 real(4),
intent(in) :: r, theta, phi
893 real(4),
dimension(3),
optional :: OCart
894 real(4),
dimension(3) :: vec, lOCart
897 if (
present(ocart)) locart = ocart
899 vec(1) = r * sin(theta) * cos(phi) - locart(1)
900 vec(2) = -r * sin(theta) * sin(phi) - locart(2)
901 vec(3) = -r * cos(theta) - locart(3)
903 END FUNCTION sphericaltocartesian4
905 FUNCTION sphericaltocartesian8 (r, theta, phi, OCart)
result (vec)
906 real(8),
intent(in) :: r, theta, phi
907 real(8),
dimension(3),
optional :: OCart
908 real(8),
dimension(3) :: vec, lOCart
911 if (
present(ocart)) locart = ocart
913 vec(1) = r * sin(theta) * cos(phi) - locart(1)
914 vec(2) = -r * sin(theta) * sin(phi) - locart(2)
915 vec(3) = -r * cos(theta) - locart(3)
917 END FUNCTION sphericaltocartesian8
919 FUNCTION sphericaltocylindrical4 (r, theta, phi, OCart)
result (vec)
920 real(4),
intent(in) :: r, theta, phi
921 real(4),
dimension(3),
optional :: OCart
922 real(4),
dimension(3) :: vec, pCart
927 pcart = sphericaltocartesian4(r, theta, phi, ocart)
928 vec = cartesiantocylindrical4(pcart(1), pcart(2), pcart(3), ocart)
930 END FUNCTION sphericaltocylindrical4
932 FUNCTION sphericaltocylindrical8 (r, theta, phi, OCart)
result (vec)
933 real(8),
intent(in) :: r, theta, phi
934 real(8),
dimension(3),
optional :: OCart
935 real(8),
dimension(3) :: vec, pCart
940 pcart = sphericaltocartesian8(r, theta, phi, ocart)
941 vec = cartesiantocylindrical8(pcart(1), pcart(2), pcart(3), ocart)
943 END FUNCTION sphericaltocylindrical8
945 FUNCTION cylindricaltospherical4 (Z, R, Phi, OCart)
result (vec)
946 real(4),
intent(in) :: Z, R, Phi
947 real(4),
dimension(3),
optional :: OCart
948 real(4),
dimension(3) :: vec, pCart
951 pcart = cylindricaltocartesian4(z, r, phi, ocart)
952 vec = cartesiantospherical4(pcart(1), pcart(2), pcart(3), ocart)
954 END FUNCTION cylindricaltospherical4
956 FUNCTION cylindricaltospherical8 (Z, R, Phi, OCart)
result (vec)
957 real(8),
intent(in) :: Z, R, Phi
958 real(8),
dimension(3),
optional :: OCart
959 real(8),
dimension(3) :: vec, pCart
962 pcart = cylindricaltocartesian8(z, r, phi, ocart)
963 vec = cartesiantospherical8(pcart(1), pcart(2), pcart(3), ocart)
965 END FUNCTION cylindricaltospherical8
968 FUNCTION currenttocartesian4 (mtype, x1, x2, x3)
result (vec)
969 integer,
intent(in) :: mtype
970 real(4),
intent(in) :: x1, x2, x3
971 real(4),
dimension(3) :: vec
973 if (mtype == mesh_types%Cartesian)
then 978 else if (mtype == mesh_types%spherical)
then 979 vec = sphericaltocartesian4(x1, x2, x3)
980 else if (mtype == mesh_types%cylindrical)
then 981 vec = cylindricaltocartesian4(x1, x2, x3)
984 END FUNCTION currenttocartesian4
986 FUNCTION currenttocartesian8 (mtype, x1, x2, x3)
result (vec)
987 integer,
intent(in) :: mtype
988 real(8),
intent(in) :: x1, x2, x3
989 real(8),
dimension(3) :: vec
991 if (mtype == mesh_types%Cartesian)
then 996 else if (mtype == mesh_types%spherical)
then 997 vec = sphericaltocartesian8(x1, x2, x3)
998 else if (mtype == mesh_types%cylindrical)
then 999 vec = cylindricaltocartesian8(x1, x2, x3)
1002 END FUNCTION currenttocartesian8
1004 FUNCTION cartesiantocurrent4 (mtype, x1, x2, x3)
result (vec)
1005 integer,
intent(in) :: mtype
1006 real(4),
intent(in) :: x1, x2, x3
1007 real(4),
dimension(3) :: vec
1009 if (mtype == mesh_types%Cartesian)
then 1014 else if (mtype == mesh_types%spherical)
then 1015 vec = cartesiantospherical4(x1, x2, x3)
1016 else if (mtype == mesh_types%cylindrical)
then 1017 vec = cartesiantocylindrical4(x1, x2, x3)
1020 END FUNCTION cartesiantocurrent4
1022 FUNCTION cartesiantocurrent8 (mtype, x1, x2, x3)
result (vec)
1023 integer,
intent(in) :: mtype
1024 real(8),
intent(in) :: x1, x2, x3
1025 real(8),
dimension(3) :: vec
1027 if (mtype == mesh_types%Cartesian)
then 1032 else if (mtype == mesh_types%spherical)
then 1033 vec = cartesiantospherical8(x1, x2, x3)
1034 else if (mtype == mesh_types%cylindrical)
then 1035 vec = cartesiantocylindrical8(x1, x2, x3)
1038 END FUNCTION cartesiantocurrent8
1044 FUNCTION volumecentroid (mesh)
result (centre)
1048 call trace_begin(
'VolumeCentroid')
1050 associate(s => mesh%s, xmn => mesh%llc_nat)
1054 centre = xmn + 0.5 * s
1056 select case (mesh%idir)
1058 centre = xmn + 0.5 * s
1062 centre = 2.0 * (3.0 * xmn**2 + 3.0 * xmn * s + s**2) &
1063 / (3.0 * (2.0 * xmn + s))
1066 select case(mesh%idir)
1070 centre = 3.0 * (4.0 * xmn**3 + 6.0 * xmn**2 * s + 4.0 * xmn * s**2 + s**3) &
1071 / (4.0 * (3.0 * xmn**2 + 3.0 * xmn * s + s**2))
1075 centre = (sin(xmn + s) - (xmn + s) * cos(xmn + s) - sin(xmn) + xmn * cos(xmn)) &
1076 / (cos(xmn) - cos(xmn + s))
1078 centre = mesh%llc_nat + 0.5 * s
1084 END FUNCTION volumecentroid
1089 FUNCTION round4(a,digit)
result (b)
1090 real(4) :: a, b, rof, q1, q2, q3
1092 real(4) :: tiny = 1.0e-32
1093 integer :: minexp = -6
1098 if (q1 <= tiny .or. digit == 0.)
then 1104 i1 = digit - int(q1 + q2 + 0.5)
1105 i2 = min(-3, digit + minexp)
1109 b = anint(a * q1 + q3) / q1
1116 FUNCTION round8(a,digit)
result (b)
1117 real(8) :: a, b, rof, q1, q2,q3
1119 real(8) :: tiny = 1.0d-99
1120 integer :: minexp = -14
1125 if (q1 <= tiny .or. digit == 0_8)
then 1130 q2 = sign(0.5_8, q1)
1131 i1 = digit - int(q1 + q2 + 0.5_8)
1132 i2 = min(-3, digit + minexp)
1136 b = anint(a * q1 + q3) / q1
Template module for mesh.