This commit is contained in:
Martin Diehl 2019-03-11 22:58:11 +01:00
parent fab813c8a1
commit 1fb1032127
3 changed files with 19 additions and 29 deletions

View File

@ -73,16 +73,16 @@ contains
function LambertCubeToBall(cube) result(ball)
use, intrinsic :: IEEE_ARITHMETIC
use prec, only: &
pInt, &
dEq0
implicit none
real(pReal), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ
real(pReal) :: T(2), c, s, q
real(pReal), dimension(2) :: T
real(pReal) :: c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal
integer(pInt), dimension(3) :: p
integer(pInt), dimension(2) :: order
integer, dimension(3) :: p
integer, dimension(2) :: order
if (maxval(abs(cube)) > AP/2.0+eps) then
ball = IEEE_value(cube,IEEE_positive_inf)
@ -135,7 +135,6 @@ pure function LambertBallToCube(xyz) result(cube)
IEEE_positive_inf, &
IEEE_value
use prec, only: &
pInt, &
dEq0
implicit none
@ -143,7 +142,7 @@ pure function LambertBallToCube(xyz) result(cube)
real(pReal), dimension(3) :: cube, xyz1, xyz3
real(pReal), dimension(2) :: Tinv, xyz2
real(pReal) :: rs, qxy, q2, sq2, q, tt
integer(pInt), dimension(3) :: p
integer, dimension(3) :: p
rs = norm2(xyz)
if (rs > R1) then
@ -192,12 +191,10 @@ end function LambertBallToCube
!> @brief determine to which pyramid a point in a cubic grid belongs
!--------------------------------------------------------------------------
pure function GetPyramidOrder(xyz)
use prec, only: &
pInt
implicit none
real(pReal),intent(in),dimension(3) :: xyz
integer(pInt), dimension(3) :: GetPyramidOrder
integer, dimension(3) :: GetPyramidOrder
if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. &
((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then

View File

@ -11,22 +11,19 @@ subroutine quit(stop_id)
use MPI, only: &
MPI_finalize
#endif
use prec, only: &
pInt
use PetscSys
use hdf5
implicit none
integer(pInt), intent(in) :: stop_id
integer, intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
integer :: hdferr
integer(pInt) :: error = 0_pInt
integer :: error
PetscErrorCode :: ierr = 0
call h5open_f(hdferr)
if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5open_f',hdferr ! prevents error if not opened yet
call h5close_f(hdferr)
if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5close_f',hdferr
call h5open_f(error)
if (error /= 0) write(6,'(a,i5)') ' Error in h5open_f ',error ! prevents error if not opened yet
call h5close_f(error)
if (error /= 0) write(6,'(a,i5)') ' Error in h5close_f ',error
call PETScFinalize(ierr)
CHKERRQ(ierr)
@ -45,8 +42,8 @@ subroutine quit(stop_id)
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0_pInt .and. ierr == 0_pInt .and. error == 0_pInt) stop 0 ! normal termination
if (stop_id == 2_pInt .and. ierr == 0_pInt .and. error == 0_pInt) stop 2 ! not all incs converged
if (stop_id == 0 .and. ierr == 0 .and. error == 0) stop 0 ! normal termination
if (stop_id == 2 .and. ierr == 0 .and. error == 0) stop 2 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit

View File

@ -383,15 +383,13 @@ end function om2eu
!> @brief convert axis angle pair to orientation matrix
!---------------------------------------------------------------------------------------------------
pure function ax2om(ax) result(om)
use prec, only: &
pInt
implicit none
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(3,3) :: om
real(pReal) :: q, c, s, omc
integer(pInt) :: i
integer :: i
c = cos(ax(4))
s = sin(ax(4))
@ -476,13 +474,12 @@ end function ax2ho
!---------------------------------------------------------------------------------------------------
pure function ho2ax(ho) result(ax)
use prec, only: &
pInt, &
dEq0
implicit none
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: ax
integer(pInt) :: i
integer :: i
real(pReal) :: hmag_squared, s, hm
real(pReal), parameter, dimension(16) :: &
tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, &
@ -519,7 +516,6 @@ end function ho2ax
!---------------------------------------------------------------------------------------------------
function om2ax(om) result(ax)
use prec, only: &
pInt, &
dEq0, &
cEq, &
dNeq0
@ -537,7 +533,7 @@ function om2ax(om) result(ax)
real(pReal), dimension(3) :: Wr, Wi
real(pReal), dimension(10) :: WORK
real(pReal), dimension(3,3) :: VR, devNull, o
integer(pInt) :: INFO, LWORK, i
integer :: INFO, LWORK, i
external :: dgeev,sgeev
@ -557,7 +553,7 @@ function om2ax(om) result(ax)
! call the eigenvalue solver
call dgeev('N','V',3,o,3,Wr,Wi,devNull,3,VR,3,WORK,LWORK,INFO)
if (INFO /= 0) call IO_error(0_pInt,ext_msg='Error in om2ax DGEEV return not zero')
if (INFO /= 0) call IO_error(0,ext_msg='Error in om2ax DGEEV return not zero')
i = maxloc(merge(1.0_pReal,0.0_pReal,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) ! poor substitute for findloc
ax(1:3) = VR(1:3,i)
where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &