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  1943998056 :         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             :