LCOV - code coverage report
Current view: | top level - atmos_phys/to_be_ccppized - linear_1d_operators.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 111 | 326 | 34.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 14 | 47 | 29.8 %

              Line data    Source code
           1             : module linear_1d_operators
           2             :
           3             :     ! This module provides the type "TriDiagOp" to represent operators on a 1D
           4             :     ! grid as tridiagonal matrices, and related types to represent boundary
           5             :     ! conditions.
           6             :     !
           7             :     ! The focus is on solving diffusion equations with a finite volume method
           8             :     ! in one dimension, but other utility operators are provided, e.g. a second
           9             :     ! order approximation to the first derivative.
          10             :     !
          11             :     ! In order to allow vectorization to occur, as well as to avoid unnecessary
          12             :     ! copying/reshaping of data in CAM, TriDiagOp actually represents a
          13             :     ! collection of independent operators that can be applied to collections of
          14             :     ! independent data; the innermost index is over independent systems (e.g.
          15             :     ! CAM columns).
          16             :     !
          17             :     ! A simple example:
          18             :     !   ! First derivative operator
          19             :     !   op = first_derivative(coords)
          20             :     !   ! Convert data to its derivative (extrapolate at boundaries).
          21             :     !   call op%apply(data)
          22             :     !
          23             :     ! With explicit boundary conditions:
          24             :     !   op = first_derivative(coords, &
          25             :     !          l_bndry=BoundaryFixedFlux(), &
          26             :     !          r_bndry=BoundaryFixedLayer(layer_distance))
          27             :     !   call op%apply(data, &
          28             :     !          l_cond=BoundaryFlux(flux, dt, thickness), &
          29             :     !          r_cond=BoundaryData(boundary))
          30             :     !
          31             :     ! Implicit solution example:
          32             :     !   ! Construct diffusion matrix.
          33             :     !   op = diffusion_operator(coords, d)
          34             :     !   call op%lmult_as_diag(-dt)
          35             :     !   call op%add_to_diag(1._r8)
          36             :     !   ! Decompose in order to invert the operation.
          37             :     !   decomp = TriDiagDecomp(op)
          38             :     !   ! Diffuse data for one time step (fixed flux boundaries).
          39             :     !   call decomp%left_div(data)
          40             :
          41             :     use shr_kind_mod, only: r8 => shr_kind_r8
          42             :     use shr_log_mod, only: errMsg => shr_log_errMsg
          43             :     use shr_sys_mod, only: shr_sys_abort
          44             :     use coords_1d, only: Coords1D
          45             :
          46             :     implicit none
          47             :     private
          48             :     save
          49             :
          50             :     ! Main type.
          51             :     public :: TriDiagOp
          52             :     public :: operator(+)
          53             :     public :: operator(-)
          54             :
          55             :     ! Decomposition used for inversion (left division).
          56             :     public :: TriDiagDecomp
          57             :
          58             :     ! Multiplies by 0.
          59             :     public :: zero_operator
          60             :
          61             :     ! Construct identity.
          62             :     public :: identity_operator
          63             :
          64             :     ! Produce a TriDiagOp that is simply a diagonal matrix.
          65             :     public :: diagonal_operator
          66             :
          67             :     ! For solving the diffusion-advection equation with implicit Euler.
          68             :     public :: diffusion_operator
          69             :     public :: advection_operator
          70             :
          71             :     ! Derivatives accurate to second order on a non-uniform grid.
          72             :     public :: first_derivative
          73             :     public :: second_derivative
          74             :
          75             :     ! Boundary condition types.
          76             :     public :: BoundaryType
          77             :     public :: BoundaryZero
          78             :     public :: BoundaryFirstOrder
          79             :     public :: BoundaryExtrapolate
          80             :     public :: BoundaryFixedLayer
          81             :     public :: BoundaryFixedFlux
          82             :
          83             :     ! Boundary data types.
          84             :     public :: BoundaryCond
          85             :     public :: BoundaryNoData
          86             :     public :: BoundaryData
          87             :     public :: BoundaryFlux
          88             :
          89             :     ! TriDiagOp represents operators that can work between nearest neighbors,
          90             :     ! with some extra logic at the boundaries. The implementation is a
          91             :     ! tridiagonal matrix plus boundary info.
          92             :     type :: TriDiagOp
          93             :        private
          94             :        ! The number of independent systems.
          95             :        integer, public :: nsys
          96             :        ! The size of the matrix (number of grid cells).
          97             :        integer, public :: ncel
          98             :        ! Super-, sub-, and regular diagonals.
          99             :        real(r8), allocatable :: spr(:,:)
         100             :        real(r8), allocatable :: sub(:,:)
         101             :        real(r8), allocatable :: diag(:,:)
         102             :        ! Buffers to hold boundary data; Details depend on the type of boundary
         103             :        ! being used.
         104             :        real(r8), allocatable :: left_bound(:)
         105             :        real(r8), allocatable :: right_bound(:)
         106             :      contains
         107             :        ! Applies the operator to a set of data.
         108             :        procedure :: apply => apply_tridiag
         109             :        ! Given the off-diagonal elements, fills in the diagonal so that the
         110             :        ! operator will have the constant function as an eigenvector with
         111             :        ! eigenvalue 0. This is used internally as a utility for construction of
         112             :        ! derivative operators.
         113             :        procedure :: deriv_diag => make_tridiag_deriv_diag
         114             :        ! Add/substract another tridiagonal from this one in-place (without
         115             :        ! creating a temporary object).
         116             :        procedure :: add => add_in_place_tridiag_ops
         117             :        procedure :: subtract => subtract_in_place_tridiag_ops
         118             :        ! Add input vector or scalar to the diagonal.
         119             :        procedure :: scalar_add_tridiag
         120             :        procedure :: diagonal_add_tridiag
         121             :        generic :: add_to_diag => scalar_add_tridiag, diagonal_add_tridiag
         122             :        ! Treat input vector (or scalar) as if it was the diagonal of an
         123             :        ! operator, and multiply this operator on the left by that value.
         124             :        procedure :: scalar_lmult_tridiag
         125             :        procedure :: diagonal_lmult_tridiag
         126             :        generic :: lmult_as_diag => &
         127             :             scalar_lmult_tridiag, diagonal_lmult_tridiag
         128             :        ! Deallocate and reset.
         129             :        procedure :: finalize => tridiag_finalize
         130             :     end type TriDiagOp
         131             :
         132             :     interface operator(+)
         133             :        module procedure add_tridiag_ops
         134             :     end interface operator(+)
         135             :
         136             :     interface operator(-)
         137             :        module procedure subtract_tridiag_ops
         138             :     end interface operator(-)
         139             :
         140             :     interface TriDiagOp
         141             :        module procedure new_TriDiagOp
         142             :     end interface TriDiagOp
         143             :
         144             :     !
         145             :     ! Boundary condition types for the operators.
         146             :     !
         147             :     ! Note that BoundaryFixedLayer and BoundaryFixedFlux are the only options
         148             :     ! supported for backwards operation (i.e. decomp%left_div). The others are
         149             :     ! meant for direct application only (e.g. to find a derivative).
         150             :     !
         151             :     ! BoundaryZero means that the operator fixes boundaries to 0.
         152             :     ! BoundaryFirstOrder means a one-sided approximation for the first
         153             :     !     derivative.
         154             :     ! BoundaryExtrapolate means that a second order approximation will be used,
         155             :     !     even at the boundaries. Boundary points do this by using their next-
         156             :     !     nearest neighbor to extrapolate.
         157             :     ! BoundaryFixedLayer means that there's an extra layer outside of the given
         158             :     !     grid, which must be specified when applying/inverting the operator.
         159             :     ! BoundaryFixedFlux is intended to provide a fixed-flux condition for
         160             :     !     typical advection/diffusion operators. It tweaks the edge condition
         161             :     !     to work on an input current rather than a value.
         162             :     !
         163             :     ! The different types were originally implemented through polymorphism, but
         164             :     ! PGI required this to be done via enum instead.
         165             :     integer, parameter :: zero_bndry = 0
         166             :     integer, parameter :: first_order_bndry = 1
         167             :     integer, parameter :: extrapolate_bndry = 2
         168             :     integer, parameter :: fixed_layer_bndry = 3
         169             :     integer, parameter :: fixed_flux_bndry = 4
         170             :
         171             :     type :: BoundaryType
         172             :        private
         173             :        integer :: bndry_type = fixed_flux_bndry
         174             :        real(r8), allocatable :: edge_width(:)
         175             :      contains
         176             :        procedure :: make_left
         177             :        procedure :: make_right
         178             :        procedure :: finalize => boundary_type_finalize
         179             :     end type BoundaryType
         180             :
         181             :     abstract interface
         182             :        subroutine deriv_seed(del_minus, del_plus, sub, spr)
         183             :          import :: r8
         184             :          real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:)
         185             :          real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:)
         186             :          real(r8), USE_CONTIGUOUS intent(out) :: sub(:)
         187             :          real(r8), USE_CONTIGUOUS intent(out) :: spr(:)
         188             :        end subroutine deriv_seed
         189             :     end interface
         190             :
         191             :     interface BoundaryZero
         192             :        module procedure new_BoundaryZero
         193             :     end interface BoundaryZero
         194             :
         195             :     interface BoundaryFirstOrder
         196             :        module procedure new_BoundaryFirstOrder
         197             :     end interface BoundaryFirstOrder
         198             :
         199             :     interface BoundaryExtrapolate
         200             :        module procedure new_BoundaryExtrapolate
         201             :     end interface BoundaryExtrapolate
         202             :
         203             :     interface BoundaryFixedLayer
         204             :        module procedure new_BoundaryFixedLayer
         205             :     end interface BoundaryFixedLayer
         206             :
         207             :     interface BoundaryFixedFlux
         208             :        module procedure new_BoundaryFixedFlux
         209             :     end interface BoundaryFixedFlux
         210             :
         211             :     !
         212             :     ! Data for boundary conditions themselves.
         213             :     !
         214             :     ! "No data" conditions perform extrapolation, if BoundaryExtrapolate was
         215             :     ! the boundary type used to construct the operator.
         216             :     !
         217             :     ! "Data" conditions contain extra data, which effectively extends the
         218             :     ! system with an extra cell.
         219             :     !
         220             :     ! "Flux" conditions contain prescribed fluxes.
         221             :     !
         222             :     ! The condition you can use depends on the boundary type from above that
         223             :     ! was used in the operator's construction. For BoundaryFixedLayer use
         224             :     ! BoundaryData. For BoundaryFixedFlux use BoundaryFlux. For everything
         225             :     ! else, use BoundaryNoData.
         226             :
         227             :     ! The switches using this enumeration used to be unnecessary due to use of
         228             :     ! polymorphism, but this had to be backed off due to insufficient PGI
         229             :     ! support for type extension.
         230             :     integer, parameter :: no_data_cond = 0
         231             :     integer, parameter :: data_cond = 1
         232             :     integer, parameter :: flux_cond = 2
         233             :
         234             :     type :: BoundaryCond
         235             :        private
         236             :        integer :: cond_type = no_data_cond
         237             :        real(r8), allocatable :: edge_data(:)
         238             :      contains
         239             :        procedure :: apply_left
         240             :        procedure :: apply_right
         241             :        procedure :: finalize => boundary_cond_finalize
         242             :     end type BoundaryCond
         243             :
         244             :     ! Constructors for different types of BoundaryCond.
         245             :     interface BoundaryNoData
         246             :        module procedure new_BoundaryNoData
         247             :     end interface BoundaryNoData
         248             :
         249             :     interface BoundaryData
         250             :        module procedure new_BoundaryData
         251             :     end interface BoundaryData
         252             :
         253             :     interface BoundaryFlux
         254             :        module procedure new_BoundaryFlux
         255             :     end interface BoundaryFlux
         256             :
         257             :     ! Opaque type to hold a tridiagonal matrix decomposition.
         258             :     !
         259             :     ! Method used is similar to Richtmyer and Morton (1967,pp 198-201), but
         260             :     ! the order of iteration is reversed, leading to A and C being swapped, and
         261             :     ! some differences in the indexing.
         262             :     type :: TriDiagDecomp
         263             :        private
         264             :        integer :: nsys = 0
         265             :        integer :: ncel = 0
         266             :        ! These correspond to A_k, E_k, and 1 / (B_k - A_k * E_{k+1})
         267             :        real(r8), allocatable :: ca(:,:)
         268             :        real(r8), allocatable :: ze(:,:)
         269             :        real(r8), allocatable :: dnom(:,:)
         270             :     contains
         271             :       procedure :: left_div => decomp_left_div
         272             :       procedure :: finalize => decomp_finalize
         273             :     end type TriDiagDecomp
         274             :
         275             :     interface TriDiagDecomp
         276             :        module procedure new_TriDiagDecomp
         277             :     end interface TriDiagDecomp
         278             :
         279             :     contains
         280             :
         281             :     ! Operator that sets to 0.
         282           0 :     function zero_operator(nsys, ncel) result(op)
         283             :       ! Sizes for operator.
         284             :       integer, intent(in) :: nsys, ncel
         285             :
         286             :       type(TriDiagOp) :: op
         287             :
         288           0 :       op = TriDiagOp(nsys, ncel)
         289             :
         290           0 :       op%spr = 0._r8
         291           0 :       op%sub = 0._r8
         292           0 :       op%diag = 0._r8
         293           0 :       op%left_bound = 0._r8
         294           0 :       op%right_bound = 0._r8
         295             :
         296           0 :     end function zero_operator
         297             :
         298             :     ! Operator that does nothing.
         299           0 :     function identity_operator(nsys, ncel) result(op)
         300             :       ! Sizes for operator.
         301             :       integer, intent(in) :: nsys, ncel
         302             :
         303             :       type(TriDiagOp) :: op
         304             :
         305           0 :       op = TriDiagOp(nsys, ncel)
         306             :
         307           0 :       op%spr = 0._r8
         308           0 :       op%sub = 0._r8
         309           0 :       op%diag = 1._r8
         310           0 :       op%left_bound = 0._r8
         311           0 :       op%right_bound = 0._r8
         312             :
         313           0 :     end function identity_operator
         314             :
         315             :     ! Create an operator that just does an element-wise product by some data.
         316     5956704 :     function diagonal_operator(diag) result(op)
         317             :       ! Data to multiply by.
         318             :       real(r8), USE_CONTIGUOUS intent(in) :: diag(:,:)
         319             :
         320             :       type(TriDiagOp) :: op
         321             :
         322     2978352 :       op = TriDiagOp(size(diag, 1), size(diag, 2))
         323             :
         324  1246267152 :       op%spr = 0._r8
         325  1246267152 :       op%sub = 0._r8
         326  1298977056 :       op%diag = diag
         327    49731552 :       op%left_bound = 0._r8
         328    49731552 :       op%right_bound = 0._r8
         329             :
         330     2978352 :     end function diagonal_operator
         331             :
         332             :     ! Diffusion matrix operator constructor. Given grid coordinates, a set of
         333             :     ! diffusion coefficients, and boundaries, creates a matrix corresponding
         334             :     ! to a finite volume representation of the operation:
         335             :     !
         336             :     ! d/dx (d_coef * d/dx)
         337             :     !
         338             :     ! This differs from what you would get from combining the first and second
         339             :     ! derivative operations, which would be more appropriate for a finite
         340             :     ! difference scheme that does not use grid cell averages.
         341    14891760 :     function diffusion_operator(coords, d_coef, l_bndry, r_bndry) &
         342             :          result(op)
         343             :       ! Grid cell locations.
         344             :       type(Coords1D), intent(in) :: coords
         345             :       ! Diffusion coefficient defined on interfaces.
         346             :       real(r8), USE_CONTIGUOUS intent(in) :: d_coef(:,:)
         347             :       ! Objects representing the kind of boundary on each side.
         348             :       class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry
         349             :       ! Output operator.
         350             :       type(TriDiagOp) :: op
         351             :
         352             :       ! Selectors to implement default boundary.
         353             :       class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc
         354             :       ! Fixed flux is default, no allocation/deallocation needed.
         355     7445880 :       type(BoundaryType), target :: bndry_default
         356             :
         357             :       ! Level index.
         358             :       integer :: k
         359             :
         360     7445880 :       if (present(l_bndry)) then
         361     1489176 :          l_bndry_loc => l_bndry
         362             :       else
         363     5956704 :          l_bndry_loc => bndry_default
         364             :       end if
         365             :
         366     7445880 :       if (present(r_bndry)) then
         367           0 :          r_bndry_loc => r_bndry
         368             :       else
         369     7445880 :          r_bndry_loc => bndry_default
         370             :       end if
         371             :
         372             :       ! Allocate the operator.
         373     7445880 :       op = TriDiagOp(coords%n, coords%d)
         374             :
         375             :       ! d_coef over the distance to the next cell gives you the matrix term for
         376             :       ! flux of material between cells. Dividing by cell thickness translates
         377             :       ! this to a tendency on the concentration. Hence the basic pattern is
         378             :       ! d_coef*rdst*rdel.
         379             :       !
         380             :       ! Boundary conditions for a fixed layer simply extend this by calculating
         381             :       ! the distance to the midpoint of the extra edge layer.
         382             :
         383     8935056 :       select case (l_bndry_loc%bndry_type)
         384             :       case (fixed_layer_bndry)
         385           0 :          op%left_bound = 2._r8*d_coef(:,1)*coords%rdel(:,1) / &
         386    26354952 :               (l_bndry_loc%edge_width+coords%del(:,1))
         387             :       case default
         388   100952280 :          op%left_bound = 0._r8
         389             :       end select
         390             :
         391   193592880 :       do k = 1, coords%d-1
         392  3108222000 :          op%spr(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k)
         393  3115667880 :          op%sub(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k+1)
         394             :       end do
         395             :
         396     7445880 :       select case (r_bndry_loc%bndry_type)
         397             :       case (fixed_layer_bndry)
         398           0 :          op%right_bound = 2._r8*d_coef(:,coords%d+1)*coords%rdel(:,coords%d) / &
         399           0 :               (r_bndry_loc%edge_width+coords%del(:,coords%d))
         400             :       case default
         401   124328880 :          op%right_bound = 0._r8
         402             :       end select
         403             :
         404             :       ! Above, we found all off-diagonals. Now get the diagonal.
         405     7445880 :       call op%deriv_diag()
         406             :
         407    14891760 :     end function diffusion_operator
         408             :
         409             :     ! Advection matrix operator constructor. Similar to diffusion_operator, it
         410             :     ! constructs an operator A corresponding to:
         411             :     !
         412             :     ! A y = d/dx (-v_coef * y)
         413             :     !
         414             :     ! Again, this is targeted at representing this operator acting on grid-cell
         415             :     ! averages in a finite volume scheme, rather than a literal representation.
         416           0 :     function advection_operator(coords, v_coef, l_bndry, r_bndry) &
         417             :          result(op)
         418             :       ! Grid cell locations.
         419             :       type(Coords1D), intent(in) :: coords
         420             :       ! Advection coefficient (effective velocity).
         421             :       real(r8), USE_CONTIGUOUS intent(in) :: v_coef(:,:)
         422             :       ! Objects representing the kind of boundary on each side.
         423             :       class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry
         424             :       ! Output operator.
         425             :       type(TriDiagOp) :: op
         426             :
         427             :       ! Selectors to implement default boundary.
         428             :       class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc
         429             :       ! Fixed flux is default, no allocation/deallocation needed.
         430           0 :       type(BoundaryType), target :: bndry_default
         431             :
         432             :       ! Negative derivative of v.
         433           0 :       real(r8) :: v_deriv(coords%n,coords%d)
         434             :
         435           0 :       if (present(l_bndry)) then
         436           0 :          l_bndry_loc => l_bndry
         437             :       else
         438           0 :          l_bndry_loc => bndry_default
         439             :       end if
         440             :
         441           0 :       if (present(r_bndry)) then
         442           0 :          r_bndry_loc => r_bndry
         443             :       else
         444           0 :          r_bndry_loc => bndry_default
         445             :       end if
         446             :
         447             :       ! Allocate the operator.
         448           0 :       op = TriDiagOp(coords%n, coords%d)
         449             :
         450             :       ! Construct the operator in two stages using the product rule. First
         451             :       ! create (-v * d/dx), then -dv/dx, and add the two.
         452             :       !
         453             :       ! For the first part, we want to interpolate to interfaces (weighted
         454             :       ! average involving del/2*dst), multiply by -v to get flux, then divide
         455             :       ! by cell thickness, which gives a concentration tendency:
         456             :       !
         457             :       !     (del/(2*dst))*(-v_coef)/del
         458             :       !
         459             :       ! Simplifying gives -v_coef*rdst*0.5, as seen below.
         460             :
         461           0 :       select case (l_bndry_loc%bndry_type)
         462             :       case (fixed_layer_bndry)
         463           0 :          op%left_bound = v_coef(:,1) / &
         464           0 :               (l_bndry_loc%edge_width+coords%del(:,1))
         465             :       case default
         466           0 :          op%left_bound = 0._r8
         467             :       end select
         468             :
         469           0 :       op%sub = v_coef(:,2:coords%d)*coords%rdst*0.5_r8
         470           0 :       op%spr = -op%sub
         471             :
         472           0 :       select case (r_bndry_loc%bndry_type)
         473             :       case (fixed_layer_bndry)
         474           0 :          op%right_bound = v_coef(:,coords%d+1) / &
         475           0 :               (r_bndry_loc%edge_width+coords%del(:,coords%d))
         476             :       case default
         477           0 :          op%right_bound = 0._r8
         478             :       end select
         479             :
         480             :       ! Above, we found all off-diagonals. Now get the diagonal. This must be
         481             :       ! done at this specific point, since the other half of the operator is
         482             :       ! not "derivative-like" in the sense of yielding 0 for a constant input.
         483           0 :       call op%deriv_diag()
         484             :
         485             :       ! The second half of the operator simply involves taking a first-order
         486             :       ! derivative of v. Since v is on the interfaces, just use:
         487             :       !     (v(k+1) - v(k))*rdel(k)
         488           0 :       v_deriv(:,1) = v_coef(:,2)*coords%rdel(:,1)
         489             :
         490           0 :       select case (l_bndry_loc%bndry_type)
         491             :       case (fixed_layer_bndry)
         492           0 :          v_deriv(:,1) = v_deriv(:,1) - v_coef(:,1)*coords%rdel(:,1)
         493             :       end select
         494             :
         495           0 :       v_deriv(:,2:coords%d-1) = (v_coef(:,3:coords%d) - &
         496           0 :            v_coef(:,2:coords%d-1))*coords%rdel(:,2:coords%d-1)
         497             :
         498           0 :       v_deriv(:,coords%d) = -v_coef(:,coords%d)*coords%rdel(:,coords%d)
         499             :
         500           0 :       select case (r_bndry_loc%bndry_type)
         501             :       case (fixed_layer_bndry)
         502             :          v_deriv(:,coords%d) = v_deriv(:,coords%d) &
         503           0 :               + v_coef(:,coords%d+1)*coords%del(:,coords%d)
         504             :       end select
         505             :
         506             :       ! Combine the two pieces.
         507           0 :       op%diag = op%diag - v_deriv
         508             :
         509           0 :     end function advection_operator
         510             :
         511             :     ! Second order approximation to the first and second derivatives on a non-
         512             :     ! uniform grid.
         513             :     !
         514             :     ! Both operators are constructed with the same method, except for a "seed"
         515             :     ! function that takes local distances between points to create the
         516             :     ! off-diagonal terms.
         517           0 :     function first_derivative(grid_spacing, l_bndry, r_bndry) result(op)
         518             :       ! Distances between points.
         519             :       real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
         520             :       ! Boundary conditions.
         521             :       class(BoundaryType), intent(in), optional :: l_bndry, r_bndry
         522             :       ! Output operator.
         523             :       type(TriDiagOp) :: op
         524             :
         525             :       op = deriv_op_from_seed(grid_spacing, first_derivative_seed, &
         526           0 :            l_bndry, r_bndry)
         527             :
         528           0 :     end function first_derivative
         529             :
         530           0 :     subroutine first_derivative_seed(del_minus, del_plus, sub, spr)
         531             :       ! Distances to next and previous point.
         532             :       real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:)
         533             :       real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:)
         534             :       ! Off-diagonal matrix terms.
         535             :       real(r8), USE_CONTIGUOUS intent(out) :: sub(:)
         536             :       real(r8), USE_CONTIGUOUS intent(out) :: spr(:)
         537             :
         538           0 :       real(r8) :: del_sum(size(del_plus))
         539             :
         540           0 :       del_sum = del_plus + del_minus
         541             :
         542           0 :       sub = - del_plus / (del_minus*del_sum)
         543           0 :       spr =   del_minus / (del_plus*del_sum)
         544             :
         545           0 :     end subroutine first_derivative_seed
         546             :
         547           0 :     function second_derivative(grid_spacing, l_bndry, r_bndry) result(op)
         548             :       ! Distances between points.
         549             :       real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
         550             :       ! Boundary conditions.
         551             :       class(BoundaryType), intent(in), optional :: l_bndry, r_bndry
         552             :       ! Output operator.
         553             :       type(TriDiagOp) :: op
         554             :
         555             :       op = deriv_op_from_seed(grid_spacing, second_derivative_seed, &
         556           0 :            l_bndry, r_bndry)
         557             :
         558           0 :     end function second_derivative
         559             :
         560           0 :     subroutine second_derivative_seed(del_minus, del_plus, sub, spr)
         561             :       ! Distances to next and previous point.
         562             :       real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:)
         563             :       real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:)
         564             :       ! Off-diagonal matrix terms.
         565             :       real(r8), USE_CONTIGUOUS intent(out) :: sub(:)
         566             :       real(r8), USE_CONTIGUOUS intent(out) :: spr(:)
         567             :
         568           0 :       real(r8) :: del_sum(size(del_plus))
         569             :
         570           0 :       del_sum = del_plus + del_minus
         571             :
         572           0 :       sub = 2._r8 / (del_minus*del_sum)
         573           0 :       spr = 2._r8 / (del_plus*del_sum)
         574             :
         575           0 :     end subroutine second_derivative_seed
         576             :
         577             :     ! Brains behind the first/second derivative functions.
         578           0 :     function deriv_op_from_seed(grid_spacing, seed, l_bndry, r_bndry) result(op)
         579             :       ! Distances between points.
         580             :       real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
         581             :       ! Function to locally construct matrix elements.
         582             :       procedure(deriv_seed) :: seed
         583             :       ! Boundary conditions.
         584             :       class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry
         585             :       ! Output operator.
         586             :       type(TriDiagOp) :: op
         587             :
         588             :       ! Selectors to implement default boundary.
         589             :       class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc
         590             :       ! Fixed flux is default, no allocation/deallocation needed.
         591           0 :       type(BoundaryType), target :: bndry_default
         592             :
         593             :       integer :: k
         594             :
         595           0 :       if (present(l_bndry)) then
         596           0 :          l_bndry_loc => l_bndry
         597             :       else
         598           0 :          l_bndry_loc => bndry_default
         599             :       end if
         600             :
         601           0 :       if (present(r_bndry)) then
         602           0 :          r_bndry_loc => r_bndry
         603             :       else
         604           0 :          r_bndry_loc => bndry_default
         605             :       end if
         606             :
         607             :       ! Number of grid points is one greater than the spacing.
         608           0 :       op = TriDiagOp(size(grid_spacing, 1), size(grid_spacing, 2) + 1)
         609             :
         610             :       ! Left boundary condition.
         611             :       call l_bndry_loc%make_left(grid_spacing, seed, &
         612           0 :            op%left_bound, op%spr(:,1))
         613             :
         614           0 :       do k = 2, op%ncel-1
         615             :          call seed(grid_spacing(:,k-1), grid_spacing(:,k), &
         616           0 :               op%sub(:,k-1), op%spr(:,k))
         617             :       end do
         618             :
         619             :       ! Right boundary condition.
         620             :       call r_bndry_loc%make_right(grid_spacing, seed, &
         621           0 :            op%sub(:,op%ncel-1), op%right_bound)
         622             :
         623             :       ! Above, we found all off-diagonals. Now get the diagonal.
         624           0 :       call op%deriv_diag()
         625             :
         626           0 :     end function deriv_op_from_seed
         627             :
         628             :     ! Boundary constructors. Most simply set an internal flag, but
         629             :     ! BoundaryFixedLayer accepts an argument representing the distance to the
         630             :     ! location where the extra layer is defined.
         631             :
         632           0 :     function new_BoundaryZero() result(new_bndry)
         633             :       type(BoundaryType) :: new_bndry
         634             :
         635           0 :       new_bndry%bndry_type = zero_bndry
         636             :
         637           0 :     end function new_BoundaryZero
         638             :
         639           0 :     function new_BoundaryFirstOrder() result(new_bndry)
         640             :       type(BoundaryType) :: new_bndry
         641             :
         642           0 :       new_bndry%bndry_type = first_order_bndry
         643             :
         644           0 :     end function new_BoundaryFirstOrder
         645             :
         646           0 :     function new_BoundaryExtrapolate() result(new_bndry)
         647             :       type(BoundaryType) :: new_bndry
         648             :
         649           0 :       new_bndry%bndry_type = extrapolate_bndry
         650             :
         651           0 :     end function new_BoundaryExtrapolate
         652             :
         653     2978352 :     function new_BoundaryFixedLayer(width) result(new_bndry)
         654             :       real(r8), USE_CONTIGUOUS intent(in) :: width(:)
         655             :       type(BoundaryType) :: new_bndry
         656             :
         657     1489176 :       new_bndry%bndry_type = fixed_layer_bndry
         658    24865776 :       new_bndry%edge_width = width
         659             :
         660     1489176 :     end function new_BoundaryFixedLayer
         661             :
         662           0 :     function new_BoundaryFixedFlux() result(new_bndry)
         663             :       type(BoundaryType) :: new_bndry
         664             :
         665             :       new_bndry%bndry_type = fixed_flux_bndry
         666             :
         667           0 :     end function new_BoundaryFixedFlux
         668             :
         669             :     ! The make_left and make_right methods implement the boundary conditions
         670             :     ! using an input seed.
         671             :
         672           0 :     subroutine make_left(self, grid_spacing, seed, term1, term2)
         673             :       class(BoundaryType), intent(in) :: self
         674             :       real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
         675             :       procedure(deriv_seed) :: seed
         676             :       real(r8), USE_CONTIGUOUS intent(out) :: term1(:)
         677             :       real(r8), USE_CONTIGUOUS intent(out) :: term2(:)
         678             :
         679           0 :       real(r8) :: del_plus(size(term1)), del_minus(size(term1))
         680             :
         681           0 :       select case (self%bndry_type)
         682             :       case (zero_bndry)
         683           0 :          term1 = 0._r8
         684           0 :          term2 = 0._r8
         685             :       case (first_order_bndry)
         686             :          ! To calculate to first order, just use a really huge del_minus (i.e.
         687             :          ! pretend that there's a point so far away it doesn't matter).
         688           0 :          del_plus = grid_spacing(:,1)
         689           0 :          del_minus = del_plus * 4._r8 / epsilon(1._r8)
         690           0 :          call seed(del_minus, del_plus, term1, term2)
         691             :       case (extrapolate_bndry)
         692             :          ! To extrapolate from the boundary, use distance from the nearest
         693             :          ! neighbor (as usual) and the second nearest neighbor (with a negative
         694             :          ! sign, since we are using two points on the same side).
         695           0 :          del_plus = grid_spacing(:,1)
         696           0 :          del_minus = - (grid_spacing(:,1) + grid_spacing(:,2))
         697           0 :          call seed(del_minus, del_plus, term1, term2)
         698             :       case (fixed_layer_bndry)
         699             :          ! Use edge value to extend the grid.
         700           0 :          del_plus = grid_spacing(:,1)
         701           0 :          del_minus = self%edge_width
         702           0 :          call seed(del_minus, del_plus, term1, term2)
         703             :       case (fixed_flux_bndry)
         704             :          ! Treat grid as uniform, but then zero out the contribution from data
         705             :          ! on one side (since it will be prescribed).
         706           0 :          del_plus = grid_spacing(:,1)
         707           0 :          del_minus = del_plus
         708           0 :          call seed(del_minus, del_plus, term1, term2)
         709           0 :          term1 = 0._r8
         710             :       case default
         711             :          call shr_sys_abort("Invalid boundary type at "// &
         712           0 :               errMsg(__FILE__, __LINE__))
         713             :       end select
         714             :
         715           0 :     end subroutine make_left
         716             :
         717           0 :     subroutine make_right(self, grid_spacing, seed, term1, term2)
         718             :       class(BoundaryType), intent(in) :: self
         719             :       real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
         720             :       procedure(deriv_seed) :: seed
         721             :       real(r8), USE_CONTIGUOUS intent(out) :: term1(:)
         722             :       real(r8), USE_CONTIGUOUS intent(out) :: term2(:)
         723             :
         724           0 :       real(r8) :: del_plus(size(term1)), del_minus(size(term1))
         725             :
         726           0 :       select case (self%bndry_type)
         727             :       case (zero_bndry)
         728           0 :          term1 = 0._r8
         729           0 :          term2 = 0._r8
         730             :       case (first_order_bndry)
         731             :          ! Use huge del_plus, analogous to how left boundary works.
         732           0 :          del_minus = grid_spacing(:,size(grid_spacing, 2))
         733           0 :          del_plus = del_minus * 4._r8 / epsilon(1._r8)
         734           0 :          call seed(del_minus, del_plus, term1, term2)
         735             :       case (extrapolate_bndry)
         736             :          ! Same strategy as left boundary, but reversed.
         737           0 :          del_plus = - (grid_spacing(:,size(grid_spacing, 2) - 1) + &
         738           0 :               grid_spacing(:,size(grid_spacing, 2)))
         739           0 :          del_minus = grid_spacing(:,size(grid_spacing, 2))
         740           0 :          call seed(del_minus, del_plus, term1, term2)
         741             :       case (fixed_layer_bndry)
         742             :          ! Use edge value to extend the grid.
         743           0 :          del_plus = self%edge_width
         744           0 :          del_minus = grid_spacing(:,size(grid_spacing, 2))
         745           0 :          call seed(del_minus, del_plus, term1, term2)
         746             :       case (fixed_flux_bndry)
         747             :          ! Uniform grid, but with edge zeroed.
         748           0 :          del_plus = grid_spacing(:,size(grid_spacing, 2))
         749           0 :          del_minus = del_plus
         750           0 :          call seed(del_minus, del_plus, term1, term2)
         751           0 :          term2 = 0._r8
         752             :       case default
         753             :          call shr_sys_abort("Invalid boundary type at "// &
         754           0 :               errMsg(__FILE__, __LINE__))
         755             :       end select
         756             :
         757           0 :     end subroutine make_right
         758             :
         759           0 :     subroutine boundary_type_finalize(self)
         760             :       class(BoundaryType), intent(inout) :: self
         761             :
         762           0 :       self%bndry_type = fixed_flux_bndry
         763           0 :       if (allocated(self%edge_width)) deallocate(self%edge_width)
         764             :
         765           0 :     end subroutine boundary_type_finalize
         766             :
         767             :     ! Constructor for TriDiagOp; this just sets the size and allocates
         768             :     ! arrays.
         769    10424232 :     type(TriDiagOp) function new_TriDiagOp(nsys, ncel)
         770             :
         771             :       integer, intent(in) :: nsys, ncel
         772             :
         773    10424232 :       new_TriDiagOp%nsys = nsys
         774    10424232 :       new_TriDiagOp%ncel = ncel
         775             :
         776             :       allocate(new_TriDiagOp%spr(nsys,ncel-1), &
         777             :            new_TriDiagOp%sub(nsys,ncel-1), &
         778             :            new_TriDiagOp%diag(nsys,ncel), &
         779             :            new_TriDiagOp%left_bound(nsys), &
         780   125090784 :            new_TriDiagOp%right_bound(nsys))
         781             :
         782    10424232 :     end function new_TriDiagOp
         783             :
         784             :     ! Deallocator for TriDiagOp.
         785    14891760 :     subroutine tridiag_finalize(self)
         786             :       class(TriDiagOp), intent(inout) :: self
         787             :
         788    14891760 :       self%nsys = 0
         789    14891760 :       self%ncel = 0
         790             :
         791    14891760 :       if (allocated(self%spr)) deallocate(self%spr)
         792    14891760 :       if (allocated(self%sub)) deallocate(self%sub)
         793    14891760 :       if (allocated(self%diag)) deallocate(self%diag)
         794    14891760 :       if (allocated(self%left_bound)) deallocate(self%left_bound)
         795    14891760 :       if (allocated(self%right_bound)) deallocate(self%right_bound)
         796             :
         797    14891760 :     end subroutine tridiag_finalize
         798             :
         799             :     ! Boundary condition constructors.
         800             :
         801           0 :     function new_BoundaryNoData() result(new_cond)
         802             :       type(BoundaryCond) :: new_cond
         803             :
         804             :       new_cond%cond_type = no_data_cond
         805             :       ! No edge data, so leave it unallocated.
         806             :
         807           0 :     end function new_BoundaryNoData
         808             :
         809     2978352 :     function new_BoundaryData(data) result(new_cond)
         810             :       real(r8), USE_CONTIGUOUS intent(in) :: data(:)
         811             :       type(BoundaryCond) :: new_cond
         812             :
         813     1489176 :       new_cond%cond_type = data_cond
         814    24865776 :       new_cond%edge_data = data
         815             :
         816     1489176 :     end function new_BoundaryData
         817             :
         818           0 :     function new_BoundaryFlux(flux, dt, spacing) result(new_cond)
         819             :       real(r8), USE_CONTIGUOUS intent(in) :: flux(:)
         820             :       real(r8), intent(in) :: dt
         821             :       real(r8), USE_CONTIGUOUS intent(in) :: spacing(:)
         822             :       type(BoundaryCond) :: new_cond
         823             :
         824           0 :       new_cond%cond_type = flux_cond
         825           0 :       new_cond%edge_data = flux*dt/spacing
         826             :
         827           0 :     end function new_BoundaryFlux
         828             :
         829             :     ! Application of input data.
         830             :     !
         831             :     ! When no data is input, assume that any bound term is applied to the
         832             :     ! third element in from the edge for extrapolation. Boundary conditions
         833             :     ! that don't need any edge data at all can then simply set the boundary
         834             :     ! terms to 0.
         835             :
         836     1489176 :     function apply_left(self, bound_term, array) result(delta_edge)
         837             :       class(BoundaryCond), intent(in) :: self
         838             :       real(r8), USE_CONTIGUOUS intent(in) :: bound_term(:)
         839             :       real(r8), USE_CONTIGUOUS intent(in) :: array(:,:)
         840             :       real(r8) :: delta_edge(size(array, 1))
         841             :
         842     1489176 :       select case (self%cond_type)
         843             :       case (no_data_cond)
         844           0 :          delta_edge = bound_term*array(:,3)
         845             :       case (data_cond)
         846    24865776 :          delta_edge = bound_term*self%edge_data
         847             :       case (flux_cond)
         848           0 :          delta_edge = self%edge_data
         849             :       case default
         850             :          call shr_sys_abort("Invalid boundary condition at "// &
         851     1489176 :               errMsg(__FILE__, __LINE__))
         852             :       end select
         853             :
         854     1489176 :     end function apply_left
         855             :
         856           0 :     function apply_right(self, bound_term, array) result(delta_edge)
         857             :       class(BoundaryCond), intent(in) :: self
         858             :       real(r8), USE_CONTIGUOUS intent(in) :: bound_term(:)
         859             :       real(r8), USE_CONTIGUOUS intent(in) :: array(:,:)
         860             :       real(r8) :: delta_edge(size(array, 1))
         861             :
         862           0 :       select case (self%cond_type)
         863             :       case (no_data_cond)
         864           0 :          delta_edge = bound_term*array(:,size(array, 2)-2)
         865             :       case (data_cond)
         866           0 :          delta_edge = bound_term*self%edge_data
         867             :       case (flux_cond)
         868           0 :          delta_edge = self%edge_data
         869             :       case default
         870             :          call shr_sys_abort("Invalid boundary condition at "// &
         871           0 :               errMsg(__FILE__, __LINE__))
         872             :       end select
         873             :
         874           0 :     end function apply_right
         875             :
         876           0 :     subroutine boundary_cond_finalize(self)
         877             :       class(BoundaryCond), intent(inout) :: self
         878             :
         879           0 :       self%cond_type = no_data_cond
         880           0 :       if (allocated(self%edge_data)) deallocate(self%edge_data)
         881             :
         882           0 :     end subroutine boundary_cond_finalize
         883             :
         884             :     ! Apply an operator and return the new data.
         885           0 :     function apply_tridiag(self, array, l_cond, r_cond) result(output)
         886             :       ! Operator to apply.
         887             :       class(TriDiagOp), intent(in) :: self
         888             :       ! Data to act on.
         889             :       real(r8), USE_CONTIGUOUS intent(in) :: array(:,:)
         890             :       ! Objects representing boundary conditions.
         891             :       class(BoundaryCond), target, intent(in), optional :: l_cond, r_cond
         892             :       ! Function result.
         893             :       real(r8) :: output(size(array, 1), size(array, 2))
         894             :
         895             :       ! Local objects to implement default.
         896             :       class(BoundaryCond), pointer :: l_cond_loc, r_cond_loc
         897             :       ! Default state is no data, no allocation/deallocation needed.
         898           0 :       type(BoundaryCond), target :: cond_default
         899             :
         900             :       ! Level index.
         901             :       integer :: k
         902             :
         903           0 :       if (present(l_cond)) then
         904           0 :          l_cond_loc => l_cond
         905             :       else
         906           0 :          l_cond_loc => cond_default
         907             :       end if
         908             :
         909           0 :       if (present(r_cond)) then
         910           0 :          r_cond_loc => r_cond
         911             :       else
         912           0 :          r_cond_loc => cond_default
         913             :       end if
         914             :
         915             :       ! Left boundary.
         916           0 :       output(:,1) = self%diag(:,1)*array(:,1) + &
         917           0 :            self%spr(:,1)*array(:,2) + &
         918           0 :            l_cond_loc%apply_left(self%left_bound, array)
         919             :
         920           0 :       do k = 2, self%ncel-1
         921           0 :          output(:,k) = &
         922           0 :               self%sub(:,k-1)*array(:,k-1) + &
         923           0 :               self%diag(:,k)*array(:,k  ) + &
         924           0 :               self%spr(:,k)*array(:,k+1)
         925             :       end do
         926             :
         927             :       ! Right boundary.
         928           0 :       output(:,self%ncel) = &
         929           0 :            self%sub(:,self%ncel-1)*array(:,self%ncel-1) + &
         930           0 :            self%diag(:,self%ncel)*array(:,self%ncel) + &
         931           0 :            r_cond_loc%apply_right(self%right_bound, array)
         932             :
         933           0 :     end function apply_tridiag
         934             :
         935             :     ! Fill in the diagonal for a TriDiagOp for a derivative operator, where
         936             :     ! the off diagonal elements are already filled in.
         937     7445880 :     subroutine make_tridiag_deriv_diag(self)
         938             :
         939             :       class(TriDiagOp), intent(inout) :: self
         940             :
         941             :       ! If a derivative operator operates on a constant function, it must
         942             :       ! return 0 everywhere. To force this, make sure that all rows add to
         943             :       ! zero in the matrix.
         944  3115667880 :       self%diag(:,:self%ncel-1) = - self%spr
         945   124328880 :       self%diag(:,self%ncel) = - self%right_bound
         946   124328880 :       self%diag(:,1) = self%diag(:,1) - self%left_bound
         947  3115667880 :       self%diag(:,2:) = self%diag(:,2:) - self%sub
         948             :
         949     7445880 :     end subroutine make_tridiag_deriv_diag
         950             :
         951             :     ! Sum two TriDiagOp objects into a new one; this is just the addition of
         952             :     ! all the entries.
         953           0 :     function add_tridiag_ops(op1, op2) result(new_op)
         954             :
         955             :       type(TriDiagOp), intent(in) :: op1, op2
         956             :       type(TriDiagOp) :: new_op
         957             :
         958           0 :       new_op = op1
         959             :
         960           0 :       call new_op%add(op2)
         961             :
         962           0 :     end function add_tridiag_ops
         963             :
         964     2978352 :     subroutine add_in_place_tridiag_ops(self, other)
         965             :
         966             :       class(TriDiagOp), intent(inout) :: self
         967             :       class(TriDiagOp), intent(in) :: other
         968             :
         969  1249245504 :       self%spr = self%spr + other%spr
         970  1249245504 :       self%sub = self%sub + other%sub
         971  1298977056 :       self%diag = self%diag + other%diag
         972             :
         973    52709904 :       self%left_bound = self%left_bound + other%left_bound
         974    52709904 :       self%right_bound = self%right_bound + other%right_bound
         975             :
         976     2978352 :     end subroutine add_in_place_tridiag_ops
         977             :
         978             :     ! Subtract two TriDiagOp objects.
         979           0 :     function subtract_tridiag_ops(op1, op2) result(new_op)
         980             :
         981             :       type(TriDiagOp), intent(in) :: op1, op2
         982             :       type(TriDiagOp) :: new_op
         983             :
         984           0 :       new_op = op1
         985             :
         986           0 :       call new_op%subtract(op2)
         987             :
         988           0 :     end function subtract_tridiag_ops
         989             :
         990             :     ! Subtract two TriDiagOp objects.
         991           0 :     subroutine subtract_in_place_tridiag_ops(self, other)
         992             :
         993             :       class(TriDiagOp), intent(inout) :: self
         994             :       class(TriDiagOp), intent(in) :: other
         995             :
         996           0 :       self%spr = self%spr - other%spr
         997           0 :       self%sub = self%sub - other%sub
         998           0 :       self%diag = self%diag - other%diag
         999             :
        1000           0 :       self%left_bound = self%left_bound - other%left_bound
        1001           0 :       self%right_bound = self%right_bound - other%right_bound
        1002             :
        1003           0 :     end subroutine subtract_in_place_tridiag_ops
        1004             :
        1005             :     ! Equivalent to adding a multiple of the identity.
        1006     7445880 :     subroutine scalar_add_tridiag(self, constant)
        1007             :
        1008             :       class(TriDiagOp), intent(inout) :: self
        1009             :       real(r8), intent(in) :: constant
        1010             :
        1011  3239996760 :       self%diag = self%diag + constant
        1012             :
        1013     7445880 :     end subroutine scalar_add_tridiag
        1014             :
        1015             :     ! Equivalent to adding the diagonal operator constructed from diag_array.
        1016           0 :     subroutine diagonal_add_tridiag(self, diag_array)
        1017             :
        1018             :       class(TriDiagOp), intent(inout) :: self
        1019             :       real(r8), USE_CONTIGUOUS intent(in) :: diag_array(:,:)
        1020             :
        1021           0 :       self%diag = self%diag + diag_array
        1022             :
        1023           0 :     end subroutine diagonal_add_tridiag
        1024             :
        1025             :     ! Multiply a scalar by an array.
        1026     7445880 :     subroutine scalar_lmult_tridiag(self, constant)
        1027             :
        1028             :       class(TriDiagOp), intent(inout) :: self
        1029             :       real(r8), intent(in) :: constant
        1030             :
        1031  3115667880 :       self%spr = self%spr * constant
        1032  3115667880 :       self%sub = self%sub * constant
        1033  3239996760 :       self%diag = self%diag * constant
        1034             :
        1035   124328880 :       self%left_bound = self%left_bound * constant
        1036   124328880 :       self%right_bound = self%right_bound * constant
        1037             :
        1038     7445880 :     end subroutine scalar_lmult_tridiag
        1039             :
        1040             :     ! Multiply in an array as if it contained the entries of a diagonal matrix
        1041             :     ! being multiplied from the left.
        1042           0 :     subroutine diagonal_lmult_tridiag(self, diag_array)
        1043             :
        1044             :       class(TriDiagOp), intent(inout) :: self
        1045             :       real(r8), USE_CONTIGUOUS intent(in) :: diag_array(:,:)
        1046             :
        1047           0 :       self%spr = self%spr * diag_array(:,:self%ncel-1)
        1048           0 :       self%sub = self%sub * diag_array(:,2:)
        1049           0 :       self%diag = self%diag * diag_array(:,:)
        1050             :
        1051           0 :       self%left_bound = self%left_bound * diag_array(:,1)
        1052           0 :       self%right_bound = self%right_bound * diag_array(:,self%ncel)
        1053             :
        1054           0 :     end subroutine diagonal_lmult_tridiag
        1055             :
        1056             :     ! Decomposition constructor
        1057             :     !
        1058             :     ! The equation to be solved later (with left_div) is:
        1059             :     !     - A(k)*q(k+1) + B(k)*q(k) - C(k)*q(k-1) = D(k)
        1060             :     !
        1061             :     ! The solution (effectively via LU decomposition) has the form:
        1062             :     !     E(k) = C(k) / (B(k) - A(k)*E(k+1))
        1063             :     !     F(k) = (D(k) + A(k)*F(k+1)) / (B(k) - A(k)*E(k+1))
        1064             :     !     q(k) = E(k) * q(k-1) + F(k)
        1065             :     !
        1066             :     ! Unlike Richtmyer and Morton, E and F are defined by iterating backward
        1067             :     ! down to level 1, and then q iterates forward.
        1068             :     !
        1069             :     ! E can be calculated and stored now, without knowing D.
        1070             :     ! To calculate F later, we store A and the denominator.
        1071     7445880 :     function new_TriDiagDecomp(op, graft_decomp) result(decomp)
        1072             :       type(TriDiagOp), intent(in) :: op
        1073             :       type(TriDiagDecomp), intent(in), optional :: graft_decomp
        1074             :
        1075             :       type(TriDiagDecomp) :: decomp
        1076             :
        1077             :       integer :: k
        1078             :
        1079     7445880 :       if (present(graft_decomp)) then
        1080           0 :          decomp%nsys = graft_decomp%nsys
        1081           0 :          decomp%ncel = graft_decomp%ncel
        1082             :       else
        1083     7445880 :          decomp%nsys = op%nsys
        1084     7445880 :          decomp%ncel = op%ncel
        1085             :       end if
        1086             :
        1087             :       ! Simple allocation with no error checking.
        1088    29783520 :       allocate(decomp%ca(decomp%nsys,decomp%ncel))
        1089    22337640 :       allocate(decomp%dnom(decomp%nsys,decomp%ncel))
        1090    22337640 :       allocate(decomp%ze(decomp%nsys,decomp%ncel))
        1091             :
        1092             :       ! decomp%ca is simply the negative of the tridiagonal's superdiagonal.
        1093  3115667880 :       decomp%ca(:,:op%ncel-1) = -op%spr
        1094   124328880 :       decomp%ca(:,op%ncel) = -op%right_bound
        1095             :
        1096     7445880 :       if (present(graft_decomp)) then
        1097             :          ! Copy in graft_decomp beyond op%ncel.
        1098           0 :          decomp%ca(:,op%ncel+1:) = graft_decomp%ca(:,op%ncel+1:)
        1099           0 :          decomp%dnom(:,op%ncel+1:) = graft_decomp%dnom(:,op%ncel+1:)
        1100           0 :          decomp%ze(:,op%ncel+1:) = graft_decomp%ze(:,op%ncel+1:)
        1101             :          ! Fill in dnom edge value.
        1102           0 :          decomp%dnom(:,op%ncel) = 1._r8 / (op%diag(:,op%ncel) - &
        1103           0 :               decomp%ca(:,op%ncel)*decomp%ze(:,op%ncel+1))
        1104             :       else
        1105             :          ! If no grafting, the edge value of dnom comes from the diagonal.
        1106   124328880 :          decomp%dnom(:,op%ncel) = 1._r8 / op%diag(:,op%ncel)
        1107             :       end if
        1108             :
        1109   193592880 :       do k = op%ncel - 1, 1, -1
        1110  3108222000 :          decomp%ze(:,k+1)   = - op%sub(:,k) * decomp%dnom(:,k+1)
        1111             :          decomp%dnom(:,k) = 1._r8 / &
        1112  3115667880 :               (op%diag(:,k) - decomp%ca(:,k)*decomp%ze(:,k+1))
        1113             :       end do
        1114             :
        1115             :       ! Don't multiply edge level by denom, because we want to leave it up to
        1116             :       ! the BoundaryCond object to decide what this means in left_div.
        1117   131774760 :       decomp%ze(:,1) = -op%left_bound
        1118             :
        1119     7445880 :     end function new_TriDiagDecomp
        1120             :
        1121             :     ! Left-division (multiplication by inverse) using a decomposed operator.
        1122             :     !
        1123             :     ! See the comment above for the constructor for a quick explanation of the
        1124             :     ! intermediate variables. The "q" argument is "D(k)" on input and "q(k)" on
        1125             :     ! output.
        1126    14891760 :     subroutine decomp_left_div(decomp, q, l_cond, r_cond)
        1127             :
        1128             :       ! Decomposed matrix.
        1129             :       class(TriDiagDecomp), intent(in) :: decomp
        1130             :       ! Data to left-divide by the matrix.
        1131             :       real(r8), USE_CONTIGUOUS intent(inout) :: q(:,:)
        1132             :       ! Objects representing boundary conditions.
        1133             :       class(BoundaryCond), intent(in), optional :: l_cond, r_cond
        1134             :
        1135             :       ! "F" from the equation above.
        1136    29783520 :       real(r8) :: zf(decomp%nsys,decomp%ncel)
        1137             :
        1138             :       ! Level index.
        1139             :       integer :: k
        1140             :
        1141             :       ! Include boundary conditions.
        1142    14891760 :       if (present(l_cond)) then
        1143    26354952 :          q(:,1) = q(:,1) + l_cond%apply_left(decomp%ze(:,1), q)
        1144             :       end if
        1145             :
        1146    14891760 :       if (present(r_cond)) then
        1147           0 :          q(:,decomp%ncel) = q(:,decomp%ncel) + &
        1148           0 :               r_cond%apply_right(decomp%ca(:,decomp%ncel), q)
        1149             :       end if
        1150             :
        1151   248657760 :       zf(:,decomp%ncel) = q(:,decomp%ncel) * decomp%dnom(:,decomp%ncel)
        1152             :
        1153   387185760 :       do k = decomp%ncel - 1, 1, -1
        1154  6231335760 :          zf(:,k) = (q(:,k) + decomp%ca(:,k)*zf(:,k+1)) * decomp%dnom(:,k)
        1155             :       end do
        1156             :
        1157             :       ! Perform back substitution
        1158             :
        1159   263549520 :       q(:,1) = zf(:,1)
        1160             :
        1161   387185760 :       do k = 2, decomp%ncel
        1162  6231335760 :          q(:,k) = zf(:,k) + decomp%ze(:,k)*q(:,k-1)
        1163             :       end do
        1164             :
        1165    14891760 :     end subroutine decomp_left_div
        1166             :
        1167             :     ! Decomposition deallocation.
        1168     7445880 :     subroutine decomp_finalize(decomp)
        1169             :       class(TriDiagDecomp), intent(inout) :: decomp
        1170             :
        1171     7445880 :       decomp%nsys = 0
        1172     7445880 :       decomp%ncel = 0
        1173             :
        1174     7445880 :       if (allocated(decomp%ca)) deallocate(decomp%ca)
        1175     7445880 :       if (allocated(decomp%dnom)) deallocate(decomp%dnom)
        1176     7445880 :       if (allocated(decomp%ze)) deallocate(decomp%ze)
        1177             :
        1178     7445880 :     end subroutine decomp_finalize
        1179             :
        1180           0 : end module linear_1d_operators
        1181             :

Generated by: LCOV version 1.14