LCOV - code coverage report
Current view: | top level - atmos_phys/to_be_ccppized - vertical_diffusion_solver.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 21 | 25 | 84.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 1 | 100.0 %

              Line data    Source code
           1             : module vertical_diffusion_solver
           2             :
           3             :     implicit none
           4             :     private
           5             :
           6             :     public :: fin_vol_solve
           7             :
           8             :     contains
           9             :
          10             :     ! Designed to solve the equation:
          11             :     !
          12             :     ! w * dq/dt = d/dp (D q' - v q) + c q
          13             :     !
          14             :     ! where q is a grid-cell average, and p is the vertical coordinate
          15             :     ! (presumably pressure).
          16             :     !
          17             :     ! In this function, coef_q_weight == w, coef_q_diff == D,
          18             :     ! coef_q_adv == v, and coef_q == c. All these are optional; omitting a
          19             :     ! coefficient is equivalent to setting the entire array to 0.
          20             :     !
          21             :     ! coef_q_diff and coef_q_adv are defined at the level interfaces, while
          22             :     ! coef_q and coef_q_weight are grid-cell averages.
          23             :
          24     4467528 :     function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_adv, &
          25     4467528 :         coef_q_weight, upper_bndry, lower_bndry, l_cond, r_cond)  result(solution)
          26             :
          27             :         use linear_1d_operators, only: &
          28             :         zero_operator,               &
          29             :         diagonal_operator,           &
          30             :         diffusion_operator,          &
          31             :         advection_operator,          &
          32             :         BoundaryType,                &
          33             :         TriDiagDecomp,               &
          34             :         TriDiagOp,                   &
          35             :         BoundaryCond,                &
          36             :         operator(+)
          37             :         use shr_kind_mod,        only: r8 => shr_kind_r8
          38             :         use coords_1d,           only: Coords1D
          39             :
          40             :         ! ---------------------- !
          41             :         ! Input-Output Arguments !
          42             :         ! ---------------------- !
          43             :
          44             :         ! Time step.
          45             :         real(r8), intent(in) :: dt
          46             :         ! Grid spacings.
          47             :         type(Coords1D), intent(in) :: p
          48             :
          49             :         integer,  intent(in)    :: ncols
          50             :         integer,  intent(in)    :: pver
          51             :         real(r8), intent(in)    :: toSolve(ncols,pver)
          52             :
          53             :         ! Coefficients for diffusion and advection.
          54             :         !
          55             :         ! The sizes must be consistent among all the coefficients that are
          56             :         ! actually present, i.e. coef_q_diff and coef_q_adv should be one level
          57             :         ! bigger than coef_q and coef_q_weight, and have the same column number.
          58             :         real(r8), contiguous, intent(in), optional :: coef_q(:,:),      &
          59             :                                                       coef_q_diff(:,:), &
          60             :                                                       coef_q_adv(:,:),  &
          61             :                                                       coef_q_weight(:,:)
          62             :
          63             :         ! Boundary conditions (optional, default to 0 flux through boundary).
          64             :         class(BoundaryType), intent(in), optional :: upper_bndry, lower_bndry
          65             :
          66             :         ! Objects representing boundary conditions.
          67             :         class(BoundaryCond), intent(in), optional :: l_cond, r_cond
          68             :
          69             :         real(r8) :: solution(ncols,pver)
          70             :
          71             :         ! decomposition.
          72     4467528 :         type(TriDiagDecomp) :: decomp
          73             :
          74             :         ! --------------- !
          75             :         ! Local Variables !
          76             :         ! --------------- !
          77             :
          78             :         ! Operator objects.
          79     4467528 :         type(TriDiagOp) :: add_term
          80     4467528 :         type(TriDiagOp) :: net_operator
          81             :
          82             :         ! ----------------------- !
          83             :         ! Main Computation Begins !
          84             :         ! ----------------------- !
          85             :
          86             :         ! A diffusion term is probably present, so start with that. Otherwise
          87             :         ! start with an operator of all 0s.
          88             :
          89     4467528 :         if (present(coef_q_diff)) then
          90     4467528 :             net_operator = diffusion_operator(p, coef_q_diff, upper_bndry, lower_bndry)
          91             :         else
          92           0 :             net_operator = zero_operator(p%n, p%d)
          93             :         end if
          94             :
          95             :         ! Constant term (damping).
          96     4467528 :         if (present(coef_q)) then
          97     2978352 :             add_term = diagonal_operator(coef_q)
          98     2978352 :             call net_operator%add(add_term)
          99             :         end if
         100             :
         101             :         ! Effective advection.
         102     4467528 :         if (present(coef_q_adv)) then
         103           0 :             add_term = advection_operator(p, coef_q_adv, upper_bndry, lower_bndry)
         104           0 :             call net_operator%add(add_term)
         105             :         end if
         106             :
         107             :         ! We want I-dt*(w^-1)*A for a single time step, implicit method, where
         108             :         ! A is the right-hand-side operator (i.e. what net_operator is now).
         109     4467528 :         if (present(coef_q_weight)) then
         110           0 :             call net_operator%lmult_as_diag(-dt/coef_q_weight)
         111             :         else
         112     4467528 :             call net_operator%lmult_as_diag(-dt)
         113             :         end if
         114     4467528 :         call net_operator%add_to_diag(1._r8)
         115             :
         116             :         ! Decompose
         117     4467528 :         decomp = TriDiagDecomp(net_operator)
         118  6942019032 :         solution = toSolve
         119             :
         120     4467528 :         call decomp%left_div(solution(:ncols, :), l_cond=l_cond, r_cond=r_cond)
         121             :
         122             :         ! Ensure local objects are deallocated.
         123     4467528 :         call decomp%finalize()
         124     4467528 :         call net_operator%finalize()
         125     4467528 :         call add_term%finalize()
         126             :
         127     8935056 :     end function fin_vol_solve
         128             :
         129             : end module vertical_diffusion_solver
         130             :

Generated by: LCOV version 1.14