pure subroutine rk2_average(ms)
!! State <- 0.5 * (state_0 + state) for h_layer, u_face_x_layer,
!! v_face_y_layer. Tracer averages are done by the caller.
type(multilayer_state_t), intent(inout) :: ms
integer :: i, j, k, nx, ny, nz, nx_face, ny_uface, nx_vface, ny_face
nx = size(ms%h_layer, 1)
ny = size(ms%h_layer, 2)
nz = ms%nz_ml
nx_face = size(ms%u_face_x_layer, 1)
ny_uface = size(ms%u_face_x_layer, 2)
nx_vface = size(ms%v_face_y_layer, 1)
ny_face = size(ms%v_face_y_layer, 2)
do concurrent(k=1:nz, j=1:ny, i=1:nx)
ms%h_layer(i, j, k) = 0.5_wp*(ms%h_layer0(i, j, k) + ms%h_layer(i, j, k))
end do
do concurrent(k=1:nz, j=1:ny_uface, i=1:nx_face)
ms%u_face_x_layer(i, j, k) = 0.5_wp*(ms%u_face_x_layer0(i, j, k) + &
ms%u_face_x_layer(i, j, k))
end do
do concurrent(k=1:nz, j=1:ny_face, i=1:nx_vface)
ms%v_face_y_layer(i, j, k) = 0.5_wp*(ms%v_face_y_layer0(i, j, k) + &
ms%v_face_y_layer(i, j, k))
end do
end subroutine rk2_average