diff --git a/src/gr/gr.f90 b/src/gr/gr.f90 index 2cb200efa..fd441c03b 100644 --- a/src/gr/gr.f90 +++ b/src/gr/gr.f90 @@ -21,7 +21,7 @@ module pure subroutine gr_getaccb_ns_body(self, system, param) real(DP) :: rmag, rdotv, vmag2 integer(I4B) :: i - associate(n => self%nbody, cb => system%cb, inv_c2 => param%inv_c2) + associate(n => self%nbody, cb => system%cb, inv_c2 => param%inv_c2) if (n == 0) return do i = 1, n rmag = norm2(self%xh(:,i)) @@ -92,10 +92,10 @@ module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh) ! Internals real(DP) :: vmag2, rmag, grterm - associate(c2 => param%inv_c2) + associate(inv_c2 => param%inv_c2) vmag2 = dot_product(pv(:), pv(:)) rmag = norm2(xh(:)) - grterm = 1.0_DP - c2 * (0.5_DP * vmag2 + 3 * mu / rmag) + grterm = 1.0_DP - inv_c2 * (0.5_DP * vmag2 + 3 * mu / rmag) vh(:) = pv(:) * grterm end associate return @@ -128,52 +128,52 @@ module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv) integer(I4B), parameter :: MAXITER = 50 real(DP),parameter :: TOL = 1.0e-12_DP - associate (c2 => param%inv_c2) + associate(inv_c2 => param%inv_c2) pv(1:NDIM) = vh(1:NDIM) ! Initial guess rterm = 3 * mu / norm2(xh(:)) v2 = dot_product(vh(:), vh(:)) do n = 1, MAXITER pv2 = dot_product(pv(:), pv(:)) - G = 1.0_DP - c2 * (0.5_DP * pv2 + rterm) + G = 1.0_DP - inv_c2 * (0.5_DP * pv2 + rterm) F(:) = pv(:) * G - vh(:) if (abs(sum(F) / v2 ) < TOL) exit ! Root found - + ! Calculate the Jacobian do k = 1, NDIM do i = 1, NDIM if (i == k) then - J(i,k) = G - c2 * pv(k) + J(i,k) = G - inv_c2 * pv(k) else - J(i,k) = - c2 * pv(k) + J(i,k) = -inv_c2 * pv(k) end if end do end do - + ! Inverse of the Jacobian det = J(1,1) * (J(3,3) * J(2,2) - J(3,2) * J(2,3)) det = det - J(2,1) * (J(3,3) * J(1,2)-J(3,2) * J(1,3)) det = det + J(3,1) * (J(2,3) * J(1,2)-J(2,2) * J(1,3)) - + Jinv(1,1) = J(3,3) * J(2,2) - J(3,2) * J(2,3) Jinv(1,2) = -(J(3,3) * J(1,2) - J(3,2) * J(1,3)) Jinv(1,3) = J(2,3) * J(1,2) - J(2,2) * J(1,3) - + Jinv(2,1) = -(J(3,3) * J(2,1) - J(3,1) * J(2,3)) Jinv(2,2) = J(3,3) * J(1,1) - J(3,1) * J(1,3) Jinv(2,3) = -(J(2,3) * J(1,1) - J(2,1) * J(1,3)) - + Jinv(3,1) = J(3,2) * J(2,1) - J(3,1) * J(2,2) Jinv(3,2) = -(J(3,2) * J(1,1) - J(3,1) * J(1,2)) Jinv(3,3) = J(2,2) * J(1,1) - J(2,1) * J(1,2) - + Jinv = Jinv * det - + do i = 1, NDIM - pv(i) = pv(i) - dot_product(Jinv(i,:) ,F(:)) + pv(i) = pv(i) - dot_product(Jinv(i,:), F(:)) end do end do - end associate + return end subroutine gr_vel2pseudovel diff --git a/src/io/io.f90 b/src/io/io.f90 index 8a6ea4ae2..dc479da5f 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -29,10 +29,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) character(STRMAX) :: line !! Line of the input file character (len=:), allocatable :: line_trim,param_name, param_value !! Strings used to parse the param file character(*),parameter :: linefmt = '(A)' !! Format code for simple text string - integer(I4B) :: integrator !! Symbolic name of integrator being used - integrator = v_list(1) ! Parse the file line by line, extracting tokens then matching them up with known parameters if possible + do read(unit = unit, fmt = linefmt, iostat = iostat, end = 1) line line_trim = trim(adjustl(line)) @@ -214,7 +213,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) write(*,*) "ENC_OUT = ",trim(adjustl(self%encounter_file)) write(*,*) "EXTRA_FORCE = ",self%lextra_force write(*,*) "BIG_DISCARD = ",self%lbig_discard - if (self%lenergy) write(*,*) "ENERGY = ",self%lenergy + write(*,*) "ENERGY = ",self%lenergy write(*,*) "MU2KG = ",self%MU2KG write(*,*) "TU2S = ",self%TU2S write(*,*) "DU2M = ",self%DU2M @@ -232,21 +231,23 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) self%inv_c2 = einsteinC * self%TU2S / self%DU2M self%inv_c2 = (self%inv_c2)**(-2) - if (integrator == RMVS) then - if (.not.self%lclose) then - write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.' - iostat = -1 - return + associate(integrator => v_list(1)) + if (integrator == RMVS) then + if (.not.self%lclose) then + write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.' + iostat = -1 + return + end if end if - end if - - ! Determine if the GR flag is set correctly for this integrator - select case(integrator) - case(WHM, RMVS) - write(*,*) "GR = ", self%lgr - case default - write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.' - end select + + ! Determine if the GR flag is set correctly for this integrator + select case(integrator) + case(WHM, RMVS) + write(*,*) "GR = ", self%lgr + case default + write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.' + end select + end associate iostat = 0 @@ -682,8 +683,8 @@ module subroutine io_read_param_in(self, param_file_name) !! Adapted from Martin Duncan's Swift routine io_init_param.f implicit none ! Arguments - class(swiftest_parameters),intent(out) :: self !! Current run configuration parameters - character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) + class(swiftest_parameters),intent(inout) :: self !! Current run configuration parameters + character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) ! Internals integer(I4B), parameter :: LUN = 7 !! Unit number of input file integer(I4B) :: ierr = 0 !! Input error code diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 0e8a0a3f5..6c81a1a11 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -488,8 +488,8 @@ end subroutine io_read_cb_in module subroutine io_read_param_in(self, param_file_name) implicit none - class(swiftest_parameters), intent(out) :: self !! Current run configuration parameters - character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) + class(swiftest_parameters), intent(inout) :: self !! Current run configuration parameters + character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) end subroutine io_read_param_in module subroutine io_read_frame_body(self, iu, param, form, ierr)