module common public :: initial,euler,output ! 2 on Unix, 8 on Mac integer, public, parameter :: double = 8 real (kind = double), public :: t,dt real (kind = double), public, dimension (2) :: pos,vel real (kind = double), public, parameter :: pi = 3.141569265358979323846 ! astronomical units real (kind = double), public, parameter :: gm = 4.0*pi*pi contains subroutine initial(nshow) integer, intent (out) :: nshow t = 0.0 print *, "time step = " read *, dt print *, "number of time steps between output = " read *, nshow print *, "initial x position = " read *, pos(1) pos(2) = 0 ! initial y-position vel(1) = 0 ! initial x-velocity print *, "initial y-velocity = " read *, vel(2) end subroutine initial subroutine euler() ! Euler-Cromer algorithm real (kind = double), dimension (2) :: accel real (kind = double) :: r2,r3 integer :: i r2 = pos(1)*pos(1) + pos(2)*pos(2) r3 = r2*sqrt(r2) do i = 1,2 accel(i) = -gm*pos(i)/r3 vel(i) = vel(i) + accel(i)*dt pos(i) = pos(i) + vel(i)*dt end do t = t + dt end subroutine euler subroutine output() print *, pos(1),pos(2) end subroutine output end module common program planet ! planetary motion use common integer :: nshow,counter call initial(nshow) call output() counter = 0 do if (t > 2) then exit end if call euler() counter = counter + 1 if (modulo(counter,nshow) == 0) then call output() ! coordinates of "earth" end if end do end program planet