Initializing water height in the tsunami simulator

Recall the initialization part of our tsunami simulator from listing 2.13, which is repeated in the following listing for your convenience, where we set the initial conditions.

Listing 3.15 Initializing the tsunami simulator

integer, parameter :: icenter = 25        ❶
real, parameter :: decay = 0.02           ❷
 
do concurrent(i = 1:grid_size)            ❸
  h(i) = exp(-decay * (i - icenter)**2)   ❸
end do                                    ❸

❶ The index at which the water perturbation will be centered

❷ The parameter that governs the width of the perturbation

❸ Sets the array values in a loop

The second part of our refactor involves defining the initialization in an external procedure so that the initial state can be set and changed more easily, by just doing

call set_gaussian(h, icenter, decay)

Now isn’t this much nicer? The algorithm is abstracted away, and if we want to change the initial parameters, we just change the values that we pass as input arguments. This subroutine will need to modify h in-place, so we’ll declare it as an intent(in out) argument, as in the “Exercise 1” sidebar. The following listing provides the complete subroutine.

Listing 3.16 A subroutine to initialize an array with a Gaussian shape

subroutine set_gaussian(x, icenter, decay)
  real, intent(in out) :: x(:)               ❶
  integer, intent(in) :: icenter             ❷
  real, intent(in) :: decay                  ❷
  integer :: i
  do concurrent(i = 1:size(x))               ❸
    x(i) = exp(-decay * (i - icenter)**2)    ❸
  end do                                     ❸
end subroutine set_gaussian

❶ One-dimensional array as input (and output) argument

❷ Input parameters for the perturbation position and shape

❸ Loops over array elements and sets their values

Like the diff function we implemented earlier, we’ll place this subroutine after the contains statement in the main program.