diff --git a/to_be_ccppized/coords_1d.F90 b/to_be_ccppized/coords_1d.F90 new file mode 100644 index 00000000..21475de0 --- /dev/null +++ b/to_be_ccppized/coords_1d.F90 @@ -0,0 +1,152 @@ +module coords_1d + + ! This module defines the Coords1D type, which is intended to to cache + ! commonly used information derived from a collection of sets of 1-D + ! coordinates. + + use shr_kind_mod, only: r8 => shr_kind_r8 + + implicit none + private + save + + public :: Coords1D + + type :: Coords1D + ! Number of sets of coordinates in the object. + integer :: n = 0 + ! Number of coordinates in each set. + integer :: d = 0 + + ! All fields below will be allocated with first dimension "n". + ! The second dimension is d+1 for ifc, d for mid, del, and rdel, and + ! d-1 for dst and rdst. + + ! Cell interface coordinates. + real(r8), allocatable :: ifc(:,:) + ! Coordinates at cell mid-points. + real(r8), allocatable :: mid(:,:) + ! Width of cells. + real(r8), allocatable :: del(:,:) + ! Distance between cell midpoints. + real(r8), allocatable :: dst(:,:) + ! Reciprocals: 1/del and 1/dst. + real(r8), allocatable :: rdel(:,:) + real(r8), allocatable :: rdst(:,:) + contains + procedure :: section + procedure :: finalize + end type Coords1D + + interface Coords1D + module procedure new_Coords1D_from_fields + module procedure new_Coords1D_from_int + end interface + + contains + + ! Constructor to create an object from existing data. + function new_Coords1D_from_fields(ifc, mid, del, dst, & + rdel, rdst) result(coords) + real(r8), USE_CONTIGUOUS intent(in) :: ifc(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: mid(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: del(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: dst(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: rdel(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: rdst(:,:) + type(Coords1D) :: coords + + coords = allocate_coords(size(ifc, 1), size(ifc, 2) - 1) + + coords%ifc = ifc + coords%mid = mid + coords%del = del + coords%dst = dst + coords%rdel = rdel + coords%rdst = rdst + + end function new_Coords1D_from_fields + + ! Constructor if you only have interface coordinates; derives all the other + ! fields. + function new_Coords1D_from_int(ifc) result(coords) + real(r8), USE_CONTIGUOUS intent(in) :: ifc(:,:) + type(Coords1D) :: coords + + coords = allocate_coords(size(ifc, 1), size(ifc, 2) - 1) + + coords%ifc = ifc + coords%mid = 0.5_r8 * (ifc(:,:coords%d)+ifc(:,2:)) + coords%del = coords%ifc(:,2:) - coords%ifc(:,:coords%d) + coords%dst = coords%mid(:,2:) - coords%mid(:,:coords%d-1) + coords%rdel = 1._r8/coords%del + coords%rdst = 1._r8/coords%dst + + end function new_Coords1D_from_int + + ! Create a new Coords1D object that is a subsection of some other object, + ! e.g. if you want only the first m coordinates, use d_bnds=[1, m]. + ! + ! Originally this used pointers, but it was found to actually be cheaper + ! in practice just to make a copy, especially since pointers can impede + ! optimization. + function section(self, n_bnds, d_bnds) + class(Coords1D), intent(in) :: self + integer, intent(in) :: n_bnds(2), d_bnds(2) + type(Coords1D) :: section + + section = allocate_coords(n_bnds(2)-n_bnds(1)+1, d_bnds(2)-d_bnds(1)+1) + + section%ifc = self%ifc(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)+1) + section%mid = self%mid(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) + section%del = self%del(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) + section%dst = self%dst(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)-1) + section%rdel = self%rdel(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) + section%rdst = self%rdst(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)-1) + + end function section + + ! Quick utility to get allocate each array with the correct size. + function allocate_coords(n, d) result(coords) + integer, intent(in) :: n, d + type(Coords1D) :: coords + + coords%n = n + coords%d = d + + allocate(coords%ifc(coords%n,coords%d+1)) + allocate(coords%mid(coords%n,coords%d)) + allocate(coords%del(coords%n,coords%d)) + allocate(coords%dst(coords%n,coords%d-1)) + allocate(coords%rdel(coords%n,coords%d)) + allocate(coords%rdst(coords%n,coords%d-1)) + + end function allocate_coords + + ! Deallocate and reset to initial state. + subroutine finalize(self) + class(Coords1D), intent(inout) :: self + + self%n = 0 + self%d = 0 + + call guarded_deallocate(self%ifc) + call guarded_deallocate(self%mid) + call guarded_deallocate(self%del) + call guarded_deallocate(self%dst) + call guarded_deallocate(self%rdel) + call guarded_deallocate(self%rdst) + + contains + + subroutine guarded_deallocate(array) + real(r8), allocatable :: array(:,:) + + if (allocated(array)) deallocate(array) + + end subroutine guarded_deallocate + + end subroutine finalize + +end module coords_1d + \ No newline at end of file diff --git a/to_be_ccppized/linear_1d_operators.F90 b/to_be_ccppized/linear_1d_operators.F90 new file mode 100644 index 00000000..901200f4 --- /dev/null +++ b/to_be_ccppized/linear_1d_operators.F90 @@ -0,0 +1,1181 @@ +module linear_1d_operators + + ! This module provides the type "TriDiagOp" to represent operators on a 1D + ! grid as tridiagonal matrices, and related types to represent boundary + ! conditions. + ! + ! The focus is on solving diffusion equations with a finite volume method + ! in one dimension, but other utility operators are provided, e.g. a second + ! order approximation to the first derivative. + ! + ! In order to allow vectorization to occur, as well as to avoid unnecessary + ! copying/reshaping of data in CAM, TriDiagOp actually represents a + ! collection of independent operators that can be applied to collections of + ! independent data; the innermost index is over independent systems (e.g. + ! CAM columns). + ! + ! A simple example: + ! ! First derivative operator + ! op = first_derivative(coords) + ! ! Convert data to its derivative (extrapolate at boundaries). + ! call op%apply(data) + ! + ! With explicit boundary conditions: + ! op = first_derivative(coords, & + ! l_bndry=BoundaryFixedFlux(), & + ! r_bndry=BoundaryFixedLayer(layer_distance)) + ! call op%apply(data, & + ! l_cond=BoundaryFlux(flux, dt, thickness), & + ! r_cond=BoundaryData(boundary)) + ! + ! Implicit solution example: + ! ! Construct diffusion matrix. + ! op = diffusion_operator(coords, d) + ! call op%lmult_as_diag(-dt) + ! call op%add_to_diag(1._r8) + ! ! Decompose in order to invert the operation. + ! decomp = TriDiagDecomp(op) + ! ! Diffuse data for one time step (fixed flux boundaries). + ! call decomp%left_div(data) + + use shr_kind_mod, only: r8 => shr_kind_r8 + use shr_log_mod, only: errMsg => shr_log_errMsg + use shr_sys_mod, only: shr_sys_abort + use coords_1d, only: Coords1D + + implicit none + private + save + + ! Main type. + public :: TriDiagOp + public :: operator(+) + public :: operator(-) + + ! Decomposition used for inversion (left division). + public :: TriDiagDecomp + + ! Multiplies by 0. + public :: zero_operator + + ! Construct identity. + public :: identity_operator + + ! Produce a TriDiagOp that is simply a diagonal matrix. + public :: diagonal_operator + + ! For solving the diffusion-advection equation with implicit Euler. + public :: diffusion_operator + public :: advection_operator + + ! Derivatives accurate to second order on a non-uniform grid. + public :: first_derivative + public :: second_derivative + + ! Boundary condition types. + public :: BoundaryType + public :: BoundaryZero + public :: BoundaryFirstOrder + public :: BoundaryExtrapolate + public :: BoundaryFixedLayer + public :: BoundaryFixedFlux + + ! Boundary data types. + public :: BoundaryCond + public :: BoundaryNoData + public :: BoundaryData + public :: BoundaryFlux + + ! TriDiagOp represents operators that can work between nearest neighbors, + ! with some extra logic at the boundaries. The implementation is a + ! tridiagonal matrix plus boundary info. + type :: TriDiagOp + private + ! The number of independent systems. + integer, public :: nsys + ! The size of the matrix (number of grid cells). + integer, public :: ncel + ! Super-, sub-, and regular diagonals. + real(r8), allocatable :: spr(:,:) + real(r8), allocatable :: sub(:,:) + real(r8), allocatable :: diag(:,:) + ! Buffers to hold boundary data; Details depend on the type of boundary + ! being used. + real(r8), allocatable :: left_bound(:) + real(r8), allocatable :: right_bound(:) + contains + ! Applies the operator to a set of data. + procedure :: apply => apply_tridiag + ! Given the off-diagonal elements, fills in the diagonal so that the + ! operator will have the constant function as an eigenvector with + ! eigenvalue 0. This is used internally as a utility for construction of + ! derivative operators. + procedure :: deriv_diag => make_tridiag_deriv_diag + ! Add/substract another tridiagonal from this one in-place (without + ! creating a temporary object). + procedure :: add => add_in_place_tridiag_ops + procedure :: subtract => subtract_in_place_tridiag_ops + ! Add input vector or scalar to the diagonal. + procedure :: scalar_add_tridiag + procedure :: diagonal_add_tridiag + generic :: add_to_diag => scalar_add_tridiag, diagonal_add_tridiag + ! Treat input vector (or scalar) as if it was the diagonal of an + ! operator, and multiply this operator on the left by that value. + procedure :: scalar_lmult_tridiag + procedure :: diagonal_lmult_tridiag + generic :: lmult_as_diag => & + scalar_lmult_tridiag, diagonal_lmult_tridiag + ! Deallocate and reset. + procedure :: finalize => tridiag_finalize + end type TriDiagOp + + interface operator(+) + module procedure add_tridiag_ops + end interface operator(+) + + interface operator(-) + module procedure subtract_tridiag_ops + end interface operator(-) + + interface TriDiagOp + module procedure new_TriDiagOp + end interface TriDiagOp + + ! + ! Boundary condition types for the operators. + ! + ! Note that BoundaryFixedLayer and BoundaryFixedFlux are the only options + ! supported for backwards operation (i.e. decomp%left_div). The others are + ! meant for direct application only (e.g. to find a derivative). + ! + ! BoundaryZero means that the operator fixes boundaries to 0. + ! BoundaryFirstOrder means a one-sided approximation for the first + ! derivative. + ! BoundaryExtrapolate means that a second order approximation will be used, + ! even at the boundaries. Boundary points do this by using their next- + ! nearest neighbor to extrapolate. + ! BoundaryFixedLayer means that there's an extra layer outside of the given + ! grid, which must be specified when applying/inverting the operator. + ! BoundaryFixedFlux is intended to provide a fixed-flux condition for + ! typical advection/diffusion operators. It tweaks the edge condition + ! to work on an input current rather than a value. + ! + ! The different types were originally implemented through polymorphism, but + ! PGI required this to be done via enum instead. + integer, parameter :: zero_bndry = 0 + integer, parameter :: first_order_bndry = 1 + integer, parameter :: extrapolate_bndry = 2 + integer, parameter :: fixed_layer_bndry = 3 + integer, parameter :: fixed_flux_bndry = 4 + + type :: BoundaryType + private + integer :: bndry_type = fixed_flux_bndry + real(r8), allocatable :: edge_width(:) + contains + procedure :: make_left + procedure :: make_right + procedure :: finalize => boundary_type_finalize + end type BoundaryType + + abstract interface + subroutine deriv_seed(del_minus, del_plus, sub, spr) + import :: r8 + real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:) + real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:) + real(r8), USE_CONTIGUOUS intent(out) :: sub(:) + real(r8), USE_CONTIGUOUS intent(out) :: spr(:) + end subroutine deriv_seed + end interface + + interface BoundaryZero + module procedure new_BoundaryZero + end interface BoundaryZero + + interface BoundaryFirstOrder + module procedure new_BoundaryFirstOrder + end interface BoundaryFirstOrder + + interface BoundaryExtrapolate + module procedure new_BoundaryExtrapolate + end interface BoundaryExtrapolate + + interface BoundaryFixedLayer + module procedure new_BoundaryFixedLayer + end interface BoundaryFixedLayer + + interface BoundaryFixedFlux + module procedure new_BoundaryFixedFlux + end interface BoundaryFixedFlux + + ! + ! Data for boundary conditions themselves. + ! + ! "No data" conditions perform extrapolation, if BoundaryExtrapolate was + ! the boundary type used to construct the operator. + ! + ! "Data" conditions contain extra data, which effectively extends the + ! system with an extra cell. + ! + ! "Flux" conditions contain prescribed fluxes. + ! + ! The condition you can use depends on the boundary type from above that + ! was used in the operator's construction. For BoundaryFixedLayer use + ! BoundaryData. For BoundaryFixedFlux use BoundaryFlux. For everything + ! else, use BoundaryNoData. + + ! The switches using this enumeration used to be unnecessary due to use of + ! polymorphism, but this had to be backed off due to insufficient PGI + ! support for type extension. + integer, parameter :: no_data_cond = 0 + integer, parameter :: data_cond = 1 + integer, parameter :: flux_cond = 2 + + type :: BoundaryCond + private + integer :: cond_type = no_data_cond + real(r8), allocatable :: edge_data(:) + contains + procedure :: apply_left + procedure :: apply_right + procedure :: finalize => boundary_cond_finalize + end type BoundaryCond + + ! Constructors for different types of BoundaryCond. + interface BoundaryNoData + module procedure new_BoundaryNoData + end interface BoundaryNoData + + interface BoundaryData + module procedure new_BoundaryData + end interface BoundaryData + + interface BoundaryFlux + module procedure new_BoundaryFlux + end interface BoundaryFlux + + ! Opaque type to hold a tridiagonal matrix decomposition. + ! + ! Method used is similar to Richtmyer and Morton (1967,pp 198-201), but + ! the order of iteration is reversed, leading to A and C being swapped, and + ! some differences in the indexing. + type :: TriDiagDecomp + private + integer :: nsys = 0 + integer :: ncel = 0 + ! These correspond to A_k, E_k, and 1 / (B_k - A_k * E_{k+1}) + real(r8), allocatable :: ca(:,:) + real(r8), allocatable :: ze(:,:) + real(r8), allocatable :: dnom(:,:) + contains + procedure :: left_div => decomp_left_div + procedure :: finalize => decomp_finalize + end type TriDiagDecomp + + interface TriDiagDecomp + module procedure new_TriDiagDecomp + end interface TriDiagDecomp + + contains + + ! Operator that sets to 0. + function zero_operator(nsys, ncel) result(op) + ! Sizes for operator. + integer, intent(in) :: nsys, ncel + + type(TriDiagOp) :: op + + op = TriDiagOp(nsys, ncel) + + op%spr = 0._r8 + op%sub = 0._r8 + op%diag = 0._r8 + op%left_bound = 0._r8 + op%right_bound = 0._r8 + + end function zero_operator + + ! Operator that does nothing. + function identity_operator(nsys, ncel) result(op) + ! Sizes for operator. + integer, intent(in) :: nsys, ncel + + type(TriDiagOp) :: op + + op = TriDiagOp(nsys, ncel) + + op%spr = 0._r8 + op%sub = 0._r8 + op%diag = 1._r8 + op%left_bound = 0._r8 + op%right_bound = 0._r8 + + end function identity_operator + + ! Create an operator that just does an element-wise product by some data. + function diagonal_operator(diag) result(op) + ! Data to multiply by. + real(r8), USE_CONTIGUOUS intent(in) :: diag(:,:) + + type(TriDiagOp) :: op + + op = TriDiagOp(size(diag, 1), size(diag, 2)) + + op%spr = 0._r8 + op%sub = 0._r8 + op%diag = diag + op%left_bound = 0._r8 + op%right_bound = 0._r8 + + end function diagonal_operator + + ! Diffusion matrix operator constructor. Given grid coordinates, a set of + ! diffusion coefficients, and boundaries, creates a matrix corresponding + ! to a finite volume representation of the operation: + ! + ! d/dx (d_coef * d/dx) + ! + ! This differs from what you would get from combining the first and second + ! derivative operations, which would be more appropriate for a finite + ! difference scheme that does not use grid cell averages. + function diffusion_operator(coords, d_coef, l_bndry, r_bndry) & + result(op) + ! Grid cell locations. + type(Coords1D), intent(in) :: coords + ! Diffusion coefficient defined on interfaces. + real(r8), USE_CONTIGUOUS intent(in) :: d_coef(:,:) + ! Objects representing the kind of boundary on each side. + class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + ! Selectors to implement default boundary. + class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc + ! Fixed flux is default, no allocation/deallocation needed. + type(BoundaryType), target :: bndry_default + + ! Level index. + integer :: k + + if (present(l_bndry)) then + l_bndry_loc => l_bndry + else + l_bndry_loc => bndry_default + end if + + if (present(r_bndry)) then + r_bndry_loc => r_bndry + else + r_bndry_loc => bndry_default + end if + + ! Allocate the operator. + op = TriDiagOp(coords%n, coords%d) + + ! d_coef over the distance to the next cell gives you the matrix term for + ! flux of material between cells. Dividing by cell thickness translates + ! this to a tendency on the concentration. Hence the basic pattern is + ! d_coef*rdst*rdel. + ! + ! Boundary conditions for a fixed layer simply extend this by calculating + ! the distance to the midpoint of the extra edge layer. + + select case (l_bndry_loc%bndry_type) + case (fixed_layer_bndry) + op%left_bound = 2._r8*d_coef(:,1)*coords%rdel(:,1) / & + (l_bndry_loc%edge_width+coords%del(:,1)) + case default + op%left_bound = 0._r8 + end select + + do k = 1, coords%d-1 + op%spr(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k) + op%sub(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k+1) + end do + + select case (r_bndry_loc%bndry_type) + case (fixed_layer_bndry) + op%right_bound = 2._r8*d_coef(:,coords%d+1)*coords%rdel(:,coords%d) / & + (r_bndry_loc%edge_width+coords%del(:,coords%d)) + case default + op%right_bound = 0._r8 + end select + + ! Above, we found all off-diagonals. Now get the diagonal. + call op%deriv_diag() + + end function diffusion_operator + + ! Advection matrix operator constructor. Similar to diffusion_operator, it + ! constructs an operator A corresponding to: + ! + ! A y = d/dx (-v_coef * y) + ! + ! Again, this is targeted at representing this operator acting on grid-cell + ! averages in a finite volume scheme, rather than a literal representation. + function advection_operator(coords, v_coef, l_bndry, r_bndry) & + result(op) + ! Grid cell locations. + type(Coords1D), intent(in) :: coords + ! Advection coefficient (effective velocity). + real(r8), USE_CONTIGUOUS intent(in) :: v_coef(:,:) + ! Objects representing the kind of boundary on each side. + class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + ! Selectors to implement default boundary. + class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc + ! Fixed flux is default, no allocation/deallocation needed. + type(BoundaryType), target :: bndry_default + + ! Negative derivative of v. + real(r8) :: v_deriv(coords%n,coords%d) + + if (present(l_bndry)) then + l_bndry_loc => l_bndry + else + l_bndry_loc => bndry_default + end if + + if (present(r_bndry)) then + r_bndry_loc => r_bndry + else + r_bndry_loc => bndry_default + end if + + ! Allocate the operator. + op = TriDiagOp(coords%n, coords%d) + + ! Construct the operator in two stages using the product rule. First + ! create (-v * d/dx), then -dv/dx, and add the two. + ! + ! For the first part, we want to interpolate to interfaces (weighted + ! average involving del/2*dst), multiply by -v to get flux, then divide + ! by cell thickness, which gives a concentration tendency: + ! + ! (del/(2*dst))*(-v_coef)/del + ! + ! Simplifying gives -v_coef*rdst*0.5, as seen below. + + select case (l_bndry_loc%bndry_type) + case (fixed_layer_bndry) + op%left_bound = v_coef(:,1) / & + (l_bndry_loc%edge_width+coords%del(:,1)) + case default + op%left_bound = 0._r8 + end select + + op%sub = v_coef(:,2:coords%d)*coords%rdst*0.5_r8 + op%spr = -op%sub + + select case (r_bndry_loc%bndry_type) + case (fixed_layer_bndry) + op%right_bound = v_coef(:,coords%d+1) / & + (r_bndry_loc%edge_width+coords%del(:,coords%d)) + case default + op%right_bound = 0._r8 + end select + + ! Above, we found all off-diagonals. Now get the diagonal. This must be + ! done at this specific point, since the other half of the operator is + ! not "derivative-like" in the sense of yielding 0 for a constant input. + call op%deriv_diag() + + ! The second half of the operator simply involves taking a first-order + ! derivative of v. Since v is on the interfaces, just use: + ! (v(k+1) - v(k))*rdel(k) + v_deriv(:,1) = v_coef(:,2)*coords%rdel(:,1) + + select case (l_bndry_loc%bndry_type) + case (fixed_layer_bndry) + v_deriv(:,1) = v_deriv(:,1) - v_coef(:,1)*coords%rdel(:,1) + end select + + v_deriv(:,2:coords%d-1) = (v_coef(:,3:coords%d) - & + v_coef(:,2:coords%d-1))*coords%rdel(:,2:coords%d-1) + + v_deriv(:,coords%d) = -v_coef(:,coords%d)*coords%rdel(:,coords%d) + + select case (r_bndry_loc%bndry_type) + case (fixed_layer_bndry) + v_deriv(:,coords%d) = v_deriv(:,coords%d) & + + v_coef(:,coords%d+1)*coords%del(:,coords%d) + end select + + ! Combine the two pieces. + op%diag = op%diag - v_deriv + + end function advection_operator + + ! Second order approximation to the first and second derivatives on a non- + ! uniform grid. + ! + ! Both operators are constructed with the same method, except for a "seed" + ! function that takes local distances between points to create the + ! off-diagonal terms. + function first_derivative(grid_spacing, l_bndry, r_bndry) result(op) + ! Distances between points. + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + ! Boundary conditions. + class(BoundaryType), intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + op = deriv_op_from_seed(grid_spacing, first_derivative_seed, & + l_bndry, r_bndry) + + end function first_derivative + + subroutine first_derivative_seed(del_minus, del_plus, sub, spr) + ! Distances to next and previous point. + real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:) + real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:) + ! Off-diagonal matrix terms. + real(r8), USE_CONTIGUOUS intent(out) :: sub(:) + real(r8), USE_CONTIGUOUS intent(out) :: spr(:) + + real(r8) :: del_sum(size(del_plus)) + + del_sum = del_plus + del_minus + + sub = - del_plus / (del_minus*del_sum) + spr = del_minus / (del_plus*del_sum) + + end subroutine first_derivative_seed + + function second_derivative(grid_spacing, l_bndry, r_bndry) result(op) + ! Distances between points. + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + ! Boundary conditions. + class(BoundaryType), intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + op = deriv_op_from_seed(grid_spacing, second_derivative_seed, & + l_bndry, r_bndry) + + end function second_derivative + + subroutine second_derivative_seed(del_minus, del_plus, sub, spr) + ! Distances to next and previous point. + real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:) + real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:) + ! Off-diagonal matrix terms. + real(r8), USE_CONTIGUOUS intent(out) :: sub(:) + real(r8), USE_CONTIGUOUS intent(out) :: spr(:) + + real(r8) :: del_sum(size(del_plus)) + + del_sum = del_plus + del_minus + + sub = 2._r8 / (del_minus*del_sum) + spr = 2._r8 / (del_plus*del_sum) + + end subroutine second_derivative_seed + + ! Brains behind the first/second derivative functions. + function deriv_op_from_seed(grid_spacing, seed, l_bndry, r_bndry) result(op) + ! Distances between points. + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + ! Function to locally construct matrix elements. + procedure(deriv_seed) :: seed + ! Boundary conditions. + class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + ! Selectors to implement default boundary. + class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc + ! Fixed flux is default, no allocation/deallocation needed. + type(BoundaryType), target :: bndry_default + + integer :: k + + if (present(l_bndry)) then + l_bndry_loc => l_bndry + else + l_bndry_loc => bndry_default + end if + + if (present(r_bndry)) then + r_bndry_loc => r_bndry + else + r_bndry_loc => bndry_default + end if + + ! Number of grid points is one greater than the spacing. + op = TriDiagOp(size(grid_spacing, 1), size(grid_spacing, 2) + 1) + + ! Left boundary condition. + call l_bndry_loc%make_left(grid_spacing, seed, & + op%left_bound, op%spr(:,1)) + + do k = 2, op%ncel-1 + call seed(grid_spacing(:,k-1), grid_spacing(:,k), & + op%sub(:,k-1), op%spr(:,k)) + end do + + ! Right boundary condition. + call r_bndry_loc%make_right(grid_spacing, seed, & + op%sub(:,op%ncel-1), op%right_bound) + + ! Above, we found all off-diagonals. Now get the diagonal. + call op%deriv_diag() + + end function deriv_op_from_seed + + ! Boundary constructors. Most simply set an internal flag, but + ! BoundaryFixedLayer accepts an argument representing the distance to the + ! location where the extra layer is defined. + + function new_BoundaryZero() result(new_bndry) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = zero_bndry + + end function new_BoundaryZero + + function new_BoundaryFirstOrder() result(new_bndry) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = first_order_bndry + + end function new_BoundaryFirstOrder + + function new_BoundaryExtrapolate() result(new_bndry) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = extrapolate_bndry + + end function new_BoundaryExtrapolate + + function new_BoundaryFixedLayer(width) result(new_bndry) + real(r8), USE_CONTIGUOUS intent(in) :: width(:) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = fixed_layer_bndry + new_bndry%edge_width = width + + end function new_BoundaryFixedLayer + + function new_BoundaryFixedFlux() result(new_bndry) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = fixed_flux_bndry + + end function new_BoundaryFixedFlux + + ! The make_left and make_right methods implement the boundary conditions + ! using an input seed. + + subroutine make_left(self, grid_spacing, seed, term1, term2) + class(BoundaryType), intent(in) :: self + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + procedure(deriv_seed) :: seed + real(r8), USE_CONTIGUOUS intent(out) :: term1(:) + real(r8), USE_CONTIGUOUS intent(out) :: term2(:) + + real(r8) :: del_plus(size(term1)), del_minus(size(term1)) + + select case (self%bndry_type) + case (zero_bndry) + term1 = 0._r8 + term2 = 0._r8 + case (first_order_bndry) + ! To calculate to first order, just use a really huge del_minus (i.e. + ! pretend that there's a point so far away it doesn't matter). + del_plus = grid_spacing(:,1) + del_minus = del_plus * 4._r8 / epsilon(1._r8) + call seed(del_minus, del_plus, term1, term2) + case (extrapolate_bndry) + ! To extrapolate from the boundary, use distance from the nearest + ! neighbor (as usual) and the second nearest neighbor (with a negative + ! sign, since we are using two points on the same side). + del_plus = grid_spacing(:,1) + del_minus = - (grid_spacing(:,1) + grid_spacing(:,2)) + call seed(del_minus, del_plus, term1, term2) + case (fixed_layer_bndry) + ! Use edge value to extend the grid. + del_plus = grid_spacing(:,1) + del_minus = self%edge_width + call seed(del_minus, del_plus, term1, term2) + case (fixed_flux_bndry) + ! Treat grid as uniform, but then zero out the contribution from data + ! on one side (since it will be prescribed). + del_plus = grid_spacing(:,1) + del_minus = del_plus + call seed(del_minus, del_plus, term1, term2) + term1 = 0._r8 + case default + call shr_sys_abort("Invalid boundary type at "// & + errMsg(__FILE__, __LINE__)) + end select + + end subroutine make_left + + subroutine make_right(self, grid_spacing, seed, term1, term2) + class(BoundaryType), intent(in) :: self + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + procedure(deriv_seed) :: seed + real(r8), USE_CONTIGUOUS intent(out) :: term1(:) + real(r8), USE_CONTIGUOUS intent(out) :: term2(:) + + real(r8) :: del_plus(size(term1)), del_minus(size(term1)) + + select case (self%bndry_type) + case (zero_bndry) + term1 = 0._r8 + term2 = 0._r8 + case (first_order_bndry) + ! Use huge del_plus, analogous to how left boundary works. + del_minus = grid_spacing(:,size(grid_spacing, 2)) + del_plus = del_minus * 4._r8 / epsilon(1._r8) + call seed(del_minus, del_plus, term1, term2) + case (extrapolate_bndry) + ! Same strategy as left boundary, but reversed. + del_plus = - (grid_spacing(:,size(grid_spacing, 2) - 1) + & + grid_spacing(:,size(grid_spacing, 2))) + del_minus = grid_spacing(:,size(grid_spacing, 2)) + call seed(del_minus, del_plus, term1, term2) + case (fixed_layer_bndry) + ! Use edge value to extend the grid. + del_plus = self%edge_width + del_minus = grid_spacing(:,size(grid_spacing, 2)) + call seed(del_minus, del_plus, term1, term2) + case (fixed_flux_bndry) + ! Uniform grid, but with edge zeroed. + del_plus = grid_spacing(:,size(grid_spacing, 2)) + del_minus = del_plus + call seed(del_minus, del_plus, term1, term2) + term2 = 0._r8 + case default + call shr_sys_abort("Invalid boundary type at "// & + errMsg(__FILE__, __LINE__)) + end select + + end subroutine make_right + + subroutine boundary_type_finalize(self) + class(BoundaryType), intent(inout) :: self + + self%bndry_type = fixed_flux_bndry + if (allocated(self%edge_width)) deallocate(self%edge_width) + + end subroutine boundary_type_finalize + + ! Constructor for TriDiagOp; this just sets the size and allocates + ! arrays. + type(TriDiagOp) function new_TriDiagOp(nsys, ncel) + + integer, intent(in) :: nsys, ncel + + new_TriDiagOp%nsys = nsys + new_TriDiagOp%ncel = ncel + + allocate(new_TriDiagOp%spr(nsys,ncel-1), & + new_TriDiagOp%sub(nsys,ncel-1), & + new_TriDiagOp%diag(nsys,ncel), & + new_TriDiagOp%left_bound(nsys), & + new_TriDiagOp%right_bound(nsys)) + + end function new_TriDiagOp + + ! Deallocator for TriDiagOp. + subroutine tridiag_finalize(self) + class(TriDiagOp), intent(inout) :: self + + self%nsys = 0 + self%ncel = 0 + + if (allocated(self%spr)) deallocate(self%spr) + if (allocated(self%sub)) deallocate(self%sub) + if (allocated(self%diag)) deallocate(self%diag) + if (allocated(self%left_bound)) deallocate(self%left_bound) + if (allocated(self%right_bound)) deallocate(self%right_bound) + + end subroutine tridiag_finalize + + ! Boundary condition constructors. + + function new_BoundaryNoData() result(new_cond) + type(BoundaryCond) :: new_cond + + new_cond%cond_type = no_data_cond + ! No edge data, so leave it unallocated. + + end function new_BoundaryNoData + + function new_BoundaryData(data) result(new_cond) + real(r8), USE_CONTIGUOUS intent(in) :: data(:) + type(BoundaryCond) :: new_cond + + new_cond%cond_type = data_cond + new_cond%edge_data = data + + end function new_BoundaryData + + function new_BoundaryFlux(flux, dt, spacing) result(new_cond) + real(r8), USE_CONTIGUOUS intent(in) :: flux(:) + real(r8), intent(in) :: dt + real(r8), USE_CONTIGUOUS intent(in) :: spacing(:) + type(BoundaryCond) :: new_cond + + new_cond%cond_type = flux_cond + new_cond%edge_data = flux*dt/spacing + + end function new_BoundaryFlux + + ! Application of input data. + ! + ! When no data is input, assume that any bound term is applied to the + ! third element in from the edge for extrapolation. Boundary conditions + ! that don't need any edge data at all can then simply set the boundary + ! terms to 0. + + function apply_left(self, bound_term, array) result(delta_edge) + class(BoundaryCond), intent(in) :: self + real(r8), USE_CONTIGUOUS intent(in) :: bound_term(:) + real(r8), USE_CONTIGUOUS intent(in) :: array(:,:) + real(r8) :: delta_edge(size(array, 1)) + + select case (self%cond_type) + case (no_data_cond) + delta_edge = bound_term*array(:,3) + case (data_cond) + delta_edge = bound_term*self%edge_data + case (flux_cond) + delta_edge = self%edge_data + case default + call shr_sys_abort("Invalid boundary condition at "// & + errMsg(__FILE__, __LINE__)) + end select + + end function apply_left + + function apply_right(self, bound_term, array) result(delta_edge) + class(BoundaryCond), intent(in) :: self + real(r8), USE_CONTIGUOUS intent(in) :: bound_term(:) + real(r8), USE_CONTIGUOUS intent(in) :: array(:,:) + real(r8) :: delta_edge(size(array, 1)) + + select case (self%cond_type) + case (no_data_cond) + delta_edge = bound_term*array(:,size(array, 2)-2) + case (data_cond) + delta_edge = bound_term*self%edge_data + case (flux_cond) + delta_edge = self%edge_data + case default + call shr_sys_abort("Invalid boundary condition at "// & + errMsg(__FILE__, __LINE__)) + end select + + end function apply_right + + subroutine boundary_cond_finalize(self) + class(BoundaryCond), intent(inout) :: self + + self%cond_type = no_data_cond + if (allocated(self%edge_data)) deallocate(self%edge_data) + + end subroutine boundary_cond_finalize + + ! Apply an operator and return the new data. + function apply_tridiag(self, array, l_cond, r_cond) result(output) + ! Operator to apply. + class(TriDiagOp), intent(in) :: self + ! Data to act on. + real(r8), USE_CONTIGUOUS intent(in) :: array(:,:) + ! Objects representing boundary conditions. + class(BoundaryCond), target, intent(in), optional :: l_cond, r_cond + ! Function result. + real(r8) :: output(size(array, 1), size(array, 2)) + + ! Local objects to implement default. + class(BoundaryCond), pointer :: l_cond_loc, r_cond_loc + ! Default state is no data, no allocation/deallocation needed. + type(BoundaryCond), target :: cond_default + + ! Level index. + integer :: k + + if (present(l_cond)) then + l_cond_loc => l_cond + else + l_cond_loc => cond_default + end if + + if (present(r_cond)) then + r_cond_loc => r_cond + else + r_cond_loc => cond_default + end if + + ! Left boundary. + output(:,1) = self%diag(:,1)*array(:,1) + & + self%spr(:,1)*array(:,2) + & + l_cond_loc%apply_left(self%left_bound, array) + + do k = 2, self%ncel-1 + output(:,k) = & + self%sub(:,k-1)*array(:,k-1) + & + self%diag(:,k)*array(:,k ) + & + self%spr(:,k)*array(:,k+1) + end do + + ! Right boundary. + output(:,self%ncel) = & + self%sub(:,self%ncel-1)*array(:,self%ncel-1) + & + self%diag(:,self%ncel)*array(:,self%ncel) + & + r_cond_loc%apply_right(self%right_bound, array) + + end function apply_tridiag + + ! Fill in the diagonal for a TriDiagOp for a derivative operator, where + ! the off diagonal elements are already filled in. + subroutine make_tridiag_deriv_diag(self) + + class(TriDiagOp), intent(inout) :: self + + ! If a derivative operator operates on a constant function, it must + ! return 0 everywhere. To force this, make sure that all rows add to + ! zero in the matrix. + self%diag(:,:self%ncel-1) = - self%spr + self%diag(:,self%ncel) = - self%right_bound + self%diag(:,1) = self%diag(:,1) - self%left_bound + self%diag(:,2:) = self%diag(:,2:) - self%sub + + end subroutine make_tridiag_deriv_diag + + ! Sum two TriDiagOp objects into a new one; this is just the addition of + ! all the entries. + function add_tridiag_ops(op1, op2) result(new_op) + + type(TriDiagOp), intent(in) :: op1, op2 + type(TriDiagOp) :: new_op + + new_op = op1 + + call new_op%add(op2) + + end function add_tridiag_ops + + subroutine add_in_place_tridiag_ops(self, other) + + class(TriDiagOp), intent(inout) :: self + class(TriDiagOp), intent(in) :: other + + self%spr = self%spr + other%spr + self%sub = self%sub + other%sub + self%diag = self%diag + other%diag + + self%left_bound = self%left_bound + other%left_bound + self%right_bound = self%right_bound + other%right_bound + + end subroutine add_in_place_tridiag_ops + + ! Subtract two TriDiagOp objects. + function subtract_tridiag_ops(op1, op2) result(new_op) + + type(TriDiagOp), intent(in) :: op1, op2 + type(TriDiagOp) :: new_op + + new_op = op1 + + call new_op%subtract(op2) + + end function subtract_tridiag_ops + + ! Subtract two TriDiagOp objects. + subroutine subtract_in_place_tridiag_ops(self, other) + + class(TriDiagOp), intent(inout) :: self + class(TriDiagOp), intent(in) :: other + + self%spr = self%spr - other%spr + self%sub = self%sub - other%sub + self%diag = self%diag - other%diag + + self%left_bound = self%left_bound - other%left_bound + self%right_bound = self%right_bound - other%right_bound + + end subroutine subtract_in_place_tridiag_ops + + ! Equivalent to adding a multiple of the identity. + subroutine scalar_add_tridiag(self, constant) + + class(TriDiagOp), intent(inout) :: self + real(r8), intent(in) :: constant + + self%diag = self%diag + constant + + end subroutine scalar_add_tridiag + + ! Equivalent to adding the diagonal operator constructed from diag_array. + subroutine diagonal_add_tridiag(self, diag_array) + + class(TriDiagOp), intent(inout) :: self + real(r8), USE_CONTIGUOUS intent(in) :: diag_array(:,:) + + self%diag = self%diag + diag_array + + end subroutine diagonal_add_tridiag + + ! Multiply a scalar by an array. + subroutine scalar_lmult_tridiag(self, constant) + + class(TriDiagOp), intent(inout) :: self + real(r8), intent(in) :: constant + + self%spr = self%spr * constant + self%sub = self%sub * constant + self%diag = self%diag * constant + + self%left_bound = self%left_bound * constant + self%right_bound = self%right_bound * constant + + end subroutine scalar_lmult_tridiag + + ! Multiply in an array as if it contained the entries of a diagonal matrix + ! being multiplied from the left. + subroutine diagonal_lmult_tridiag(self, diag_array) + + class(TriDiagOp), intent(inout) :: self + real(r8), USE_CONTIGUOUS intent(in) :: diag_array(:,:) + + self%spr = self%spr * diag_array(:,:self%ncel-1) + self%sub = self%sub * diag_array(:,2:) + self%diag = self%diag * diag_array(:,:) + + self%left_bound = self%left_bound * diag_array(:,1) + self%right_bound = self%right_bound * diag_array(:,self%ncel) + + end subroutine diagonal_lmult_tridiag + + ! Decomposition constructor + ! + ! The equation to be solved later (with left_div) is: + ! - A(k)*q(k+1) + B(k)*q(k) - C(k)*q(k-1) = D(k) + ! + ! The solution (effectively via LU decomposition) has the form: + ! E(k) = C(k) / (B(k) - A(k)*E(k+1)) + ! F(k) = (D(k) + A(k)*F(k+1)) / (B(k) - A(k)*E(k+1)) + ! q(k) = E(k) * q(k-1) + F(k) + ! + ! Unlike Richtmyer and Morton, E and F are defined by iterating backward + ! down to level 1, and then q iterates forward. + ! + ! E can be calculated and stored now, without knowing D. + ! To calculate F later, we store A and the denominator. + function new_TriDiagDecomp(op, graft_decomp) result(decomp) + type(TriDiagOp), intent(in) :: op + type(TriDiagDecomp), intent(in), optional :: graft_decomp + + type(TriDiagDecomp) :: decomp + + integer :: k + + if (present(graft_decomp)) then + decomp%nsys = graft_decomp%nsys + decomp%ncel = graft_decomp%ncel + else + decomp%nsys = op%nsys + decomp%ncel = op%ncel + end if + + ! Simple allocation with no error checking. + allocate(decomp%ca(decomp%nsys,decomp%ncel)) + allocate(decomp%dnom(decomp%nsys,decomp%ncel)) + allocate(decomp%ze(decomp%nsys,decomp%ncel)) + + ! decomp%ca is simply the negative of the tridiagonal's superdiagonal. + decomp%ca(:,:op%ncel-1) = -op%spr + decomp%ca(:,op%ncel) = -op%right_bound + + if (present(graft_decomp)) then + ! Copy in graft_decomp beyond op%ncel. + decomp%ca(:,op%ncel+1:) = graft_decomp%ca(:,op%ncel+1:) + decomp%dnom(:,op%ncel+1:) = graft_decomp%dnom(:,op%ncel+1:) + decomp%ze(:,op%ncel+1:) = graft_decomp%ze(:,op%ncel+1:) + ! Fill in dnom edge value. + decomp%dnom(:,op%ncel) = 1._r8 / (op%diag(:,op%ncel) - & + decomp%ca(:,op%ncel)*decomp%ze(:,op%ncel+1)) + else + ! If no grafting, the edge value of dnom comes from the diagonal. + decomp%dnom(:,op%ncel) = 1._r8 / op%diag(:,op%ncel) + end if + + do k = op%ncel - 1, 1, -1 + decomp%ze(:,k+1) = - op%sub(:,k) * decomp%dnom(:,k+1) + decomp%dnom(:,k) = 1._r8 / & + (op%diag(:,k) - decomp%ca(:,k)*decomp%ze(:,k+1)) + end do + + ! Don't multiply edge level by denom, because we want to leave it up to + ! the BoundaryCond object to decide what this means in left_div. + decomp%ze(:,1) = -op%left_bound + + end function new_TriDiagDecomp + + ! Left-division (multiplication by inverse) using a decomposed operator. + ! + ! See the comment above for the constructor for a quick explanation of the + ! intermediate variables. The "q" argument is "D(k)" on input and "q(k)" on + ! output. + subroutine decomp_left_div(decomp, q, l_cond, r_cond) + + ! Decomposed matrix. + class(TriDiagDecomp), intent(in) :: decomp + ! Data to left-divide by the matrix. + real(r8), USE_CONTIGUOUS intent(inout) :: q(:,:) + ! Objects representing boundary conditions. + class(BoundaryCond), intent(in), optional :: l_cond, r_cond + + ! "F" from the equation above. + real(r8) :: zf(decomp%nsys,decomp%ncel) + + ! Level index. + integer :: k + + ! Include boundary conditions. + if (present(l_cond)) then + q(:,1) = q(:,1) + l_cond%apply_left(decomp%ze(:,1), q) + end if + + if (present(r_cond)) then + q(:,decomp%ncel) = q(:,decomp%ncel) + & + r_cond%apply_right(decomp%ca(:,decomp%ncel), q) + end if + + zf(:,decomp%ncel) = q(:,decomp%ncel) * decomp%dnom(:,decomp%ncel) + + do k = decomp%ncel - 1, 1, -1 + zf(:,k) = (q(:,k) + decomp%ca(:,k)*zf(:,k+1)) * decomp%dnom(:,k) + end do + + ! Perform back substitution + + q(:,1) = zf(:,1) + + do k = 2, decomp%ncel + q(:,k) = zf(:,k) + decomp%ze(:,k)*q(:,k-1) + end do + + end subroutine decomp_left_div + + ! Decomposition deallocation. + subroutine decomp_finalize(decomp) + class(TriDiagDecomp), intent(inout) :: decomp + + decomp%nsys = 0 + decomp%ncel = 0 + + if (allocated(decomp%ca)) deallocate(decomp%ca) + if (allocated(decomp%dnom)) deallocate(decomp%dnom) + if (allocated(decomp%ze)) deallocate(decomp%ze) + + end subroutine decomp_finalize + +end module linear_1d_operators + \ No newline at end of file diff --git a/to_be_ccppized/vertical_diffusion_solver.F90 b/to_be_ccppized/vertical_diffusion_solver.F90 new file mode 100644 index 00000000..c49d0c45 --- /dev/null +++ b/to_be_ccppized/vertical_diffusion_solver.F90 @@ -0,0 +1,130 @@ +module vertical_diffusion_solver + + implicit none + private + + public :: fin_vol_solve + + contains + + ! Designed to solve the equation: + ! + ! w * dq/dt = d/dp (D q' - v q) + c q + ! + ! where q is a grid-cell average, and p is the vertical coordinate + ! (presumably pressure). + ! + ! In this function, coef_q_weight == w, coef_q_diff == D, + ! coef_q_adv == v, and coef_q == c. All these are optional; omitting a + ! coefficient is equivalent to setting the entire array to 0. + ! + ! coef_q_diff and coef_q_adv are defined at the level interfaces, while + ! coef_q and coef_q_weight are grid-cell averages. + + function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_adv, & + coef_q_weight, upper_bndry, lower_bndry, l_cond, r_cond) result(solution) + + use linear_1d_operators, only: & + zero_operator, & + diagonal_operator, & + diffusion_operator, & + advection_operator, & + BoundaryType, & + TriDiagDecomp, & + TriDiagOp, & + BoundaryCond, & + operator(+) + use shr_kind_mod, only: r8 => shr_kind_r8 + use coords_1d, only: Coords1D + + ! ---------------------- ! + ! Input-Output Arguments ! + ! ---------------------- ! + + ! Time step. + real(r8), intent(in) :: dt + ! Grid spacings. + type(Coords1D), intent(in) :: p + + integer, intent(in) :: ncols + integer, intent(in) :: pver + real(r8), intent(in) :: toSolve(ncols,pver) + + ! Coefficients for diffusion and advection. + ! + ! The sizes must be consistent among all the coefficients that are + ! actually present, i.e. coef_q_diff and coef_q_adv should be one level + ! bigger than coef_q and coef_q_weight, and have the same column number. + real(r8), contiguous, intent(in), optional :: coef_q(:,:), & + coef_q_diff(:,:), & + coef_q_adv(:,:), & + coef_q_weight(:,:) + + ! Boundary conditions (optional, default to 0 flux through boundary). + class(BoundaryType), intent(in), optional :: upper_bndry, lower_bndry + + ! Objects representing boundary conditions. + class(BoundaryCond), intent(in), optional :: l_cond, r_cond + + real(r8) :: solution(ncols,pver) + + ! decomposition. + type(TriDiagDecomp) :: decomp + + ! --------------- ! + ! Local Variables ! + ! --------------- ! + + ! Operator objects. + type(TriDiagOp) :: add_term + type(TriDiagOp) :: net_operator + + ! ----------------------- ! + ! Main Computation Begins ! + ! ----------------------- ! + + ! A diffusion term is probably present, so start with that. Otherwise + ! start with an operator of all 0s. + + if (present(coef_q_diff)) then + net_operator = diffusion_operator(p, coef_q_diff, upper_bndry, lower_bndry) + else + net_operator = zero_operator(p%n, p%d) + end if + + ! Constant term (damping). + if (present(coef_q)) then + add_term = diagonal_operator(coef_q) + call net_operator%add(add_term) + end if + + ! Effective advection. + if (present(coef_q_adv)) then + add_term = advection_operator(p, coef_q_adv, upper_bndry, lower_bndry) + call net_operator%add(add_term) + end if + + ! We want I-dt*(w^-1)*A for a single time step, implicit method, where + ! A is the right-hand-side operator (i.e. what net_operator is now). + if (present(coef_q_weight)) then + call net_operator%lmult_as_diag(-dt/coef_q_weight) + else + call net_operator%lmult_as_diag(-dt) + end if + call net_operator%add_to_diag(1._r8) + + ! Decompose + decomp = TriDiagDecomp(net_operator) + solution = toSolve + + call decomp%left_div(solution(:ncols, :), l_cond=l_cond, r_cond=r_cond) + + ! Ensure local objects are deallocated. + call decomp%finalize() + call net_operator%finalize() + call add_term%finalize() + + end function fin_vol_solve + +end module vertical_diffusion_solver + \ No newline at end of file