module common public :: initial,print_table,Euler_Richardson ! 2 on Unix, 8 on Mac integer, public, parameter :: double = 2 real (kind = double), public :: y,v,a,t,vt2,dt,dt_2 real (kind = double), public, parameter :: g = 9.8 contains subroutine initial(y0,t0) real (kind = double), intent (in out) :: y0,t0 real (kind = double) :: vt t0 = -0.132 ! initial time (see Table 3.1) y0 = 0 y = y0 v = 0 ! initial velocity print *, "terminal velocity = " read *, vt vt2 = vt*vt dt = 0.001 dt_2 = 0.5*dt end subroutine initial subroutine print_table(y0,nshow) real (kind = double), intent (in out) :: y0 integer, intent (out) :: nshow real (kind = double) :: show_time,distance_fallen if (t < 0) then show_time = 0.1 ! time interval between output ! choice of dt = 0.001 convenient nshow = int(show_time/dt) print "(t8,a,t20,a)", "time","displacement" print *, "" end if distance_fallen = y0 - y print "(2f13.4)", t,distance_fallen end subroutine print_table subroutine Euler_Richardson() ! Euler-Richardson method real (kind = double) :: v2,vmid,v2mid,ymid,a,amid v2 = v*v a = -g*(1 - v2/vt2) vmid = v + a*dt_2 ! velocity at midpoint ymid = y + v*dt_2 ! position at midpoint v2mid = vmid*vmid amid = -g*(1 - v2mid/vt2) ! acceleration at midpoint v = v + amid*dt y = y + vmid*dt t = t + dt end subroutine Euler_Richardson end module common program drag ! assume drag force proportional to v^2 use common real (kind = double) :: y0 integer :: nshow,counter call initial(y0,t) call print_table(y0,nshow) iterate: do if (t >= 0) then exit iterate end if call Euler_Richardson() end do iterate call print_table(y0,nshow) ! values at t = 0 counter = 0 loop: do if (t > 0.8) then exit loop end if call Euler_Richardson() counter = counter + 1 if (modulo(counter,nshow) == 0) then call print_table(y0,nshow) end if end do loop end program drag