module m_les !! Shared physics for explicit eddy-viscosity LES models. !! !! Derivatives and data movement are orchestrated here, while pointwise !! operations are dispatched to the selected computational backend. use m_base_backend, only: base_backend_t use m_common, only: dp, DIR_X, DIR_Y, DIR_Z, DIR_C, VERT, & RDR_X2Y, RDR_X2Z, RDR_Y2X, RDR_Z2X use m_config, only: les_config_t use m_field, only: field_t use m_mesh, only: mesh_t use m_tdsops, only: dirps_t, tdsops_t implicit none private public :: les_t, smagorinsky_nut, strain_rate_magnitude, & filter_width, wall_damped_mixing_length, neutral_wall_stress, & neutral_drag_coefficient type :: les_t character(len=20) :: model = 'none' real(dp) :: smagorinsky_constant = 0.14_dp logical :: wall_damping = .false. real(dp) :: wall_damping_n = 3._dp !> Wall the Mason-Thomson damping measures from, supplied by the case !> through configure_wall_damping (the ABL case uses its kappa and z0) logical :: wall_supplied = .false. real(dp) :: von_karman_constant = 0._dp real(dp) :: roughness_length = 0._dp logical :: abl_wall_boundary_enabled = .false. real(dp) :: abl_wall_sampling_height = 0._dp !! y-vertex the wall model samples the velocity at (1 is the wall). integer :: abl_wall_sample_plane = 2 class(field_t), pointer :: nut => null() class(field_t), pointer :: mixing_length_sq => null() contains procedure :: mixing_length procedure :: nut_from_gradient procedure :: compute_nut procedure :: apply_sgs_stress procedure :: configure_wall_damping procedure :: configure_abl_wall_boundary procedure :: finalise end type les_t interface les_t module procedure les_init end interface les_t contains pure function les_init(config) result(les) type(les_config_t), intent(in) :: config type(les_t) :: les les%model = config%model les%smagorinsky_constant = config%smagorinsky_constant les%wall_damping = config%wall_damping les%wall_damping_n = config%wall_damping_n end function les_init pure real(dp) function filter_width(spacing) result(delta) !! Volume-equivalent grid-filter width, Delta=(dx*dy*dz)^(1/3). real(dp), intent(in) :: spacing(3) if (any(spacing <= 0._dp)) then delta = 0._dp else delta = product(spacing)**(1._dp/3._dp) end if end function filter_width pure real(dp) function strain_rate_magnitude(velocity_gradient) result(smag) !! |S|=sqrt(2 S_ij S_ij), where S_ij=(du_i/dx_j+du_j/dx_i)/2. real(dp), intent(in) :: velocity_gradient(3, 3) real(dp) :: sij_sq sij_sq = velocity_gradient(1, 1)**2 + & velocity_gradient(2, 2)**2 + & velocity_gradient(3, 3)**2 + & 0.5_dp*(velocity_gradient(1, 2) + & velocity_gradient(2, 1))**2 + & 0.5_dp*(velocity_gradient(1, 3) + & velocity_gradient(3, 1))**2 + & 0.5_dp*(velocity_gradient(2, 3) + & velocity_gradient(3, 2))**2 smag = sqrt(2._dp*sij_sq) end function strain_rate_magnitude pure real(dp) function wall_damped_mixing_length( & delta, wall_distance, smagorinsky_constant, von_karman_constant, & exponent, roughness_length) result(length) !! Mason-Thomson blend used by Incompact3d's ABL Smagorinsky path. real(dp), intent(in) :: delta, wall_distance, smagorinsky_constant real(dp), intent(in) :: von_karman_constant, exponent, roughness_length real(dp) :: wall_scale if (delta <= 0._dp .or. wall_distance + roughness_length <= 0._dp) then length = 0._dp return end if wall_scale = von_karman_constant*(wall_distance + roughness_length)/delta length = delta*(smagorinsky_constant**(-exponent) + & wall_scale**(-exponent))**(-1._dp/exponent) end function wall_damped_mixing_length pure real(dp) function smagorinsky_nut( & velocity_gradient, length) result(nut) real(dp), intent(in) :: velocity_gradient(3, 3) real(dp), intent(in) :: length nut = length**2*strain_rate_magnitude(velocity_gradient) end function smagorinsky_nut pure real(dp) function mixing_length( & self, spacing, wall_distance) result(length) !! Smagorinsky mixing length, optionally Mason-Thomson wall-damped. !! !! Returns zero when wall damping is requested without a wall distance, !! which in turn zeroes the eddy viscosity built from it. class(les_t), intent(in) :: self real(dp), intent(in) :: spacing(3) real(dp), optional, intent(in) :: wall_distance real(dp) :: delta delta = filter_width(spacing) length = self%smagorinsky_constant*delta if (self%wall_damping) then if (.not. present(wall_distance)) then length = 0._dp return end if length = wall_damped_mixing_length( & delta, wall_distance, & self%smagorinsky_constant, self%von_karman_constant, & self%wall_damping_n, self%roughness_length & ) end if end function mixing_length pure real(dp) function nut_from_gradient( & self, velocity_gradient, spacing, wall_distance) result(nut) class(les_t), intent(in) :: self real(dp), intent(in) :: velocity_gradient(3, 3), spacing(3) real(dp), optional, intent(in) :: wall_distance if (trim(self%model) == 'none') then nut = 0._dp return end if nut = smagorinsky_nut(velocity_gradient, & self%mixing_length(spacing, wall_distance)) end function nut_from_gradient subroutine configure_wall_damping(self, kappa, roughness_length) !! Supply the wall the Mason-Thomson damping measures from. The wall !! distance is taken from the lower y boundary, so this suits cases with a !! single wall there. class(les_t), intent(inout) :: self real(dp), intent(in) :: kappa, roughness_length if (associated(self%mixing_length_sq)) & error stop 'Configure the LES wall before applying LES.' if (kappa <= 0._dp) & error stop 'LES wall damping needs a positive von Karman constant.' if (roughness_length < 0._dp) & error stop 'LES wall damping needs a non-negative roughness length.' self%von_karman_constant = kappa self%roughness_length = roughness_length self%wall_supplied = .true. end subroutine configure_wall_damping subroutine configure_abl_wall_boundary( & self, kappa, roughness_length, sampling_height, sample_plane) class(les_t), intent(inout) :: self real(dp), intent(in) :: kappa, roughness_length, sampling_height integer, intent(in) :: sample_plane if (trim(self%model) /= 'smagorinsky') & error stop 'The neutral ABL wall model requires Smagorinsky LES.' if (.not. self%wall_damping) & error stop 'The neutral ABL wall model requires LES wall damping.' if (sample_plane < 2) & error stop 'The ABL wall model must sample above the no-slip floor.' ! ABL owns the wall properties. LES consumes the same values for its ! Mason-Thomson damping and for the wall stress. call self%configure_wall_damping(kappa, roughness_length) self%abl_wall_boundary_enabled = .true. self%abl_wall_sampling_height = sampling_height self%abl_wall_sample_plane = sample_plane end subroutine configure_abl_wall_boundary subroutine compute_nut(self, backend, mesh, u, v, w, & xdirps, ydirps, zdirps) !! Compute the Smagorinsky eddy-viscosity field in DIR_X layout. class(les_t), intent(inout) :: self class(base_backend_t), intent(inout) :: backend type(mesh_t), intent(in) :: mesh class(field_t), intent(in) :: u, v, w type(dirps_t), intent(in) :: xdirps, ydirps, zdirps class(field_t), pointer :: dudx, dudy, dudz class(field_t), pointer :: dvdx, dvdy, dvdz class(field_t), pointer :: dwdx, dwdy, dwdz if (u%dir /= DIR_X .or. v%dir /= DIR_X .or. w%dir /= DIR_X) then error stop 'LES velocity fields must use DIR_X layout.' end if if (.not. associated(self%nut)) then self%nut => backend%allocator%get_block(DIR_X, VERT) end if select case (trim(self%model)) case ('none') call self%nut%fill(0._dp) return case ('smagorinsky') case default error stop 'Unsupported LES model in compute_nut.' end select if (.not. associated(self%mixing_length_sq)) then self%mixing_length_sq => backend%allocator%get_block(DIR_X, VERT) call initialise_mixing_length(self, backend, mesh) end if call compute_velocity_gradients( & backend, u, v, w, xdirps, ydirps, zdirps, & dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz) call backend%compute_smagorinsky_nut( & self%nut, self%mixing_length_sq, & dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz) call release_velocity_gradients( & backend, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz) end subroutine compute_nut subroutine apply_sgs_stress(self, backend, mesh, du, dv, dw, u, v, w, & xdirps, ydirps, zdirps) !! Add div(2*nut*S_ij) to the momentum right-hand-side fields. class(les_t), intent(inout) :: self class(base_backend_t), intent(inout) :: backend type(mesh_t), intent(in) :: mesh class(field_t), intent(inout) :: du, dv, dw class(field_t), intent(in) :: u, v, w type(dirps_t), intent(in) :: xdirps, ydirps, zdirps class(field_t), pointer :: dudx, dudy, dudz class(field_t), pointer :: dvdx, dvdy, dvdz class(field_t), pointer :: dwdx, dwdy, dwdz real(dp) :: wall_drag_coeff if (trim(self%model) == 'none') return if (trim(self%model) /= 'smagorinsky') & error stop 'Unsupported LES model in apply_sgs_stress.' if (du%dir /= DIR_X .or. dv%dir /= DIR_X .or. dw%dir /= DIR_X) then error stop 'LES momentum RHS fields must use DIR_X layout.' end if if (u%dir /= DIR_X .or. v%dir /= DIR_X .or. w%dir /= DIR_X) then error stop 'LES velocity fields must use DIR_X layout.' end if if (.not. associated(self%nut)) & self%nut => backend%allocator%get_block(DIR_X, VERT) if (.not. associated(self%mixing_length_sq)) then self%mixing_length_sq => backend%allocator%get_block(DIR_X, VERT) call initialise_mixing_length(self, backend, mesh) end if call compute_velocity_gradients( & backend, u, v, w, xdirps, ydirps, zdirps, & dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz) call backend%compute_smagorinsky_nut( & self%nut, self%mixing_length_sq, & dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz) wall_drag_coeff = 0._dp if (self%abl_wall_boundary_enabled) then wall_drag_coeff = neutral_drag_coefficient( & self%von_karman_constant, self%roughness_length, & self%abl_wall_sampling_height) end if call add_sgs_terms(backend, du, dv, dw, self%nut, & dudx, dudy, dudz, dvdx, dvdy, dvdz, & dwdx, dwdy, dwdz, xdirps, ydirps, zdirps, & self%abl_wall_boundary_enabled, u, w, & self%abl_wall_sample_plane, wall_drag_coeff) call release_velocity_gradients( & backend, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz) end subroutine apply_sgs_stress subroutine add_sgs_terms(backend, du, dv, dw, nut, & dudx, dudy, dudz, dvdx, dvdy, dvdz, & dwdx, dwdy, dwdz, xdirps, ydirps, zdirps, & abl_wall, wall_u, wall_w, & wall_sample_plane, wall_drag_coeff) class(base_backend_t), intent(inout) :: backend class(field_t), intent(inout) :: du, dv, dw class(field_t), intent(in) :: nut class(field_t), intent(in) :: dudx, dudy, dudz class(field_t), intent(in) :: dvdx, dvdy, dvdz class(field_t), intent(in) :: dwdx, dwdy, dwdz type(dirps_t), intent(in) :: xdirps, ydirps, zdirps logical, intent(in) :: abl_wall class(field_t), intent(in) :: wall_u, wall_w integer, intent(in) :: wall_sample_plane real(dp), intent(in) :: wall_drag_coeff ! On the first plane above a no-slip floor the whole SGS stress tensor is ! replaced, as Incompact3d does in sgs_mom_conservative: tau_xy and tau_yz ! take the modelled wall stress, and the remaining four components are ! cleared so the resolved gradient across the no-slip condition cannot add ! a second stress there. call add_normal_stress(backend, du, nut, dudx, xdirps, abl_wall) call add_normal_stress(backend, dv, nut, dvdy, ydirps, abl_wall) call add_normal_stress(backend, dw, nut, dwdz, zdirps, abl_wall) call add_shear_stress(backend, du, ydirps, dv, xdirps, nut, dudy, dvdx, & abl_wall, 1, wall_u, wall_w, wall_sample_plane, & wall_drag_coeff) call add_shear_stress(backend, du, zdirps, dw, xdirps, nut, dudz, dwdx, & abl_wall, 0, wall_u, wall_w, wall_sample_plane, & wall_drag_coeff) call add_shear_stress(backend, dv, zdirps, dw, ydirps, nut, dvdz, dwdy, & abl_wall, 3, wall_u, wall_w, wall_sample_plane, & wall_drag_coeff) end subroutine add_sgs_terms pure subroutine neutral_wall_stress(u_sample, w_sample, kappa, & roughness_length, sampling_height, & tau_x, tau_z) !! Neutral rough-wall drag law. Returns the wall value of the SGS stress !! tau_xy (and tau_yz), signed so that a positive sample gives a positive !! stress and hence a momentum sink once differentiated. real(dp), intent(in) :: u_sample, w_sample, kappa real(dp), intent(in) :: roughness_length, sampling_height real(dp), intent(out) :: tau_x, tau_z real(dp) :: drag_coeff, speed drag_coeff = neutral_drag_coefficient(kappa, roughness_length, & sampling_height) speed = sqrt(u_sample**2 + w_sample**2) tau_x = drag_coeff*u_sample*speed tau_z = drag_coeff*w_sample*speed end subroutine neutral_wall_stress pure real(dp) function neutral_drag_coefficient( & kappa, roughness_length, sampling_height) result(drag_coeff) !! Neutral rough-wall drag coefficient, (kappa/ln(h/z0))**2, so that !! the wall stress is drag_coeff*u*|u| at the sampling height h. real(dp), intent(in) :: kappa, roughness_length, sampling_height drag_coeff = (kappa/log(sampling_height/roughness_length))**2 end function neutral_drag_coefficient subroutine compute_velocity_gradients( & backend, u, v, w, xdirps, ydirps, zdirps, & dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz) class(base_backend_t), intent(inout) :: backend class(field_t), intent(in) :: u, v, w type(dirps_t), intent(in) :: xdirps, ydirps, zdirps class(field_t), pointer, intent(out) :: dudx, dudy, dudz class(field_t), pointer, intent(out) :: dvdx, dvdy, dvdz class(field_t), pointer, intent(out) :: dwdx, dwdy, dwdz ! The sym flag only matters at a free-slip (Neumann) boundary, which the ! compact schemes close by mirroring the field across it. There the ! component normal to the boundary is odd (it vanishes on it) and the ! tangential ones are even (zero normal gradient). So du_i/dx_j takes ! the odd operator when i = j (u_i is normal to a j-boundary) and the ! even one when i /= j (u_i is tangential to it). For the ABL this ! selects, at the free-slip lid, the odd closure for dv/dy and the even ! one for du/dy and dw/dy; x and z are periodic, where the flag has no ! effect, and they follow the same rule for consistency. This matches ! Incompact3d's npaire choice for the same derivatives. call derivative_to_x(backend, dudx, u, xdirps, sym=.false.) call derivative_to_x(backend, dvdx, v, xdirps, sym=.true.) call derivative_to_x(backend, dwdx, w, xdirps, sym=.true.) call derivative_to_x(backend, dudy, u, ydirps, sym=.true.) call derivative_to_x(backend, dvdy, v, ydirps, sym=.false.) call derivative_to_x(backend, dwdy, w, ydirps, sym=.true.) call derivative_to_x(backend, dudz, u, zdirps, sym=.true.) call derivative_to_x(backend, dvdz, v, zdirps, sym=.true.) call derivative_to_x(backend, dwdz, w, zdirps, sym=.false.) end subroutine compute_velocity_gradients subroutine release_velocity_gradients( & backend, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz) class(base_backend_t), intent(inout) :: backend class(field_t), pointer, intent(inout) :: dudx, dudy, dudz class(field_t), pointer, intent(inout) :: dvdx, dvdy, dvdz class(field_t), pointer, intent(inout) :: dwdx, dwdy, dwdz call backend%allocator%release_block(dudx) call backend%allocator%release_block(dudy) call backend%allocator%release_block(dudz) call backend%allocator%release_block(dvdx) call backend%allocator%release_block(dvdy) call backend%allocator%release_block(dvdz) call backend%allocator%release_block(dwdx) call backend%allocator%release_block(dwdy) call backend%allocator%release_block(dwdz) end subroutine release_velocity_gradients subroutine add_normal_stress(backend, rhs, nut, gradient, direction, & stamp_wall) class(base_backend_t), intent(inout) :: backend class(field_t), intent(inout) :: rhs class(field_t), intent(in) :: nut, gradient type(dirps_t), intent(in) :: direction !! When set, clear this stress on the first plane above a no-slip floor. logical, intent(in) :: stamp_wall class(field_t), pointer :: stress stress => backend%allocator%get_block(DIR_X, VERT) call backend%compute_sgs_stress( & stress, nut, gradient, gradient, 2._dp, 0._dp) if (stamp_wall) call backend%field_set_y_plane(stress, 0._dp, 2) ! tau_ii is even across a free-slip boundary in direction i. call add_stress_derivative(backend, rhs, stress, direction, sym=.true.) call backend%allocator%release_block(stress) end subroutine add_normal_stress subroutine add_shear_stress( & backend, rhs_a, direction_a, rhs_b, & direction_b, nut, gradient_a, gradient_b, & stamp_wall, wall_component, wall_u, wall_w, wall_sample_plane, & wall_drag_coeff & ) class(base_backend_t), intent(inout) :: backend class(field_t), intent(inout) :: rhs_a, rhs_b type(dirps_t), intent(in) :: direction_a, direction_b class(field_t), intent(in) :: nut, gradient_a, gradient_b !! When set, replace the floor value of this stress with the modelled !! wall stress before differentiating, so the wall flux is carried by !! the same operator that transports it in the interior. logical, intent(in) :: stamp_wall integer, intent(in) :: wall_component, wall_sample_plane !! wall_component: 1 for tau_xy, 3 for tau_yz (drag law from the local !! sampled velocity), 0 for tau_xz, which the wall does not carry class(field_t), intent(in) :: wall_u, wall_w real(dp), intent(in) :: wall_drag_coeff class(field_t), pointer :: stress stress => backend%allocator%get_block(DIR_X, VERT) call backend%compute_sgs_stress( & stress, nut, gradient_a, gradient_b, 1._dp, 1._dp) if (stamp_wall) then ! Substitute the modelled stress at the first vertex above the no-slip ! floor, as Incompact3d does (te1/th1 at index 2 in ! sgs_mom_conservative). That is where the resolved gradient across the ! no-slip condition would otherwise produce a second, spurious stress on ! top of the modelled one. The free-slip lid carries no stress, and its ! odd closure ignores the boundary value, so it is left alone. if (wall_component == 0) then call backend%field_set_y_plane(stress, 0._dp, 2) else call backend%field_set_abl_wall_stress( & stress, wall_u, wall_w, wall_sample_plane, 2, & wall_drag_coeff, wall_component) end if end if ! tau_ij (i/=j) is odd across a free-slip boundary in directions i and j. call add_stress_derivative(backend, rhs_a, stress, direction_a, & sym=.false.) call add_stress_derivative(backend, rhs_b, stress, direction_b, & sym=.false.) call backend%allocator%release_block(stress) end subroutine add_shear_stress subroutine add_stress_derivative(backend, rhs, stress, direction, sym) class(base_backend_t), intent(inout) :: backend class(field_t), intent(inout) :: rhs class(field_t), intent(in) :: stress type(dirps_t), intent(in) :: direction logical, intent(in) :: sym class(field_t), pointer :: derivative call derivative_to_x(backend, derivative, stress, direction, sym) call backend%vecadd(1._dp, derivative, 1._dp, rhs) call backend%allocator%release_block(derivative) end subroutine add_stress_derivative subroutine derivative_to_x(backend, derivative, velocity, dirps, sym) class(base_backend_t), intent(inout) :: backend class(field_t), pointer, intent(out) :: derivative class(field_t), intent(in) :: velocity type(dirps_t), target, intent(in) :: dirps !! Parity of the differentiated field across a free-slip (Neumann) !! boundary: .true. selects the even (sym) operator, .false. the odd one. !! Irrelevant for periodic and Dirichlet boundaries. logical, intent(in) :: sym class(field_t), pointer :: velocity_dir, derivative_dir class(tdsops_t), pointer :: der1st if (sym) then der1st => dirps%der1st_sym else der1st => dirps%der1st end if derivative => backend%allocator%get_block(DIR_X) select case (dirps%dir) case (DIR_X) call backend%tds_solve(derivative, velocity, der1st) case (DIR_Y) velocity_dir => backend%allocator%get_block(DIR_Y) derivative_dir => backend%allocator%get_block(DIR_Y) call backend%reorder(velocity_dir, velocity, RDR_X2Y) call backend%tds_solve(derivative_dir, velocity_dir, der1st) call backend%reorder(derivative, derivative_dir, RDR_Y2X) call backend%allocator%release_block(velocity_dir) call backend%allocator%release_block(derivative_dir) case (DIR_Z) velocity_dir => backend%allocator%get_block(DIR_Z) derivative_dir => backend%allocator%get_block(DIR_Z) call backend%reorder(velocity_dir, velocity, RDR_X2Z) call backend%tds_solve(derivative_dir, velocity_dir, der1st) call backend%reorder(derivative, derivative_dir, RDR_Z2X) call backend%allocator%release_block(velocity_dir) call backend%allocator%release_block(derivative_dir) case default error stop 'Invalid derivative direction in LES.' end select end subroutine derivative_to_x subroutine initialise_mixing_length(self, backend, mesh) class(les_t), intent(in) :: self class(base_backend_t), intent(inout) :: backend type(mesh_t), intent(in) :: mesh real(dp), allocatable :: mixing_data(:, :, :) real(dp) :: spacing(3), wall_distance, y_lower integer :: dims(3), dims_padded(3), i, j, k ! The case supplies the wall after the solver (and so this les_t) is ! built, so this one-off setup is the first point where it can be checked. if (self%wall_damping .and. .not. self%wall_supplied) & error stop 'LES wall damping needs a case that supplies the wall & &(currently only abl).' dims = mesh%get_dims(VERT) dims_padded = backend%allocator%get_padded_dims(DIR_C) allocate (mixing_data(dims_padded(1), dims_padded(2), dims_padded(3))) mixing_data = 0._dp y_lower = 0._dp if (trim(mesh%geo%stretching(2)) == 'centred') & y_lower = -0.5_dp*mesh%geo%L(2) do k = 1, dims(3) spacing(3) = spacing_at_vertex(mesh, k, 3, dims(3)) do j = 1, dims(2) spacing(2) = spacing_at_vertex(mesh, j, 2, dims(2)) wall_distance = max(mesh%geo%vert_coords(j, 2) - y_lower, 0._dp) do i = 1, dims(1) spacing(1) = spacing_at_vertex(mesh, i, 1, dims(1)) mixing_data(i, j, k) = self%mixing_length(spacing, wall_distance)**2 end do end do end do call backend%set_field_data(self%mixing_length_sq, mixing_data) deallocate (mixing_data) end subroutine initialise_mixing_length pure real(dp) function spacing_at_vertex(mesh, index, direction, n) & result(spacing) type(mesh_t), intent(in) :: mesh integer, intent(in) :: index, direction, n if (.not. mesh%geo%stretched(direction) .or. n == 1) then spacing = mesh%geo%d(direction) else if (index < n) then spacing = abs(mesh%geo%vert_coords(index + 1, direction) - & mesh%geo%vert_coords(index, direction)) else spacing = abs(mesh%geo%vert_coords(index, direction) - & mesh%geo%vert_coords(index - 1, direction)) end if end function spacing_at_vertex subroutine finalise(self, backend) class(les_t), intent(inout) :: self class(base_backend_t), intent(inout) :: backend if (associated(self%nut)) then call backend%allocator%release_block(self%nut) nullify (self%nut) end if if (associated(self%mixing_length_sq)) then call backend%allocator%release_block(self%mixing_length_sq) nullify (self%mixing_length_sq) end if end subroutine finalise end module m_les