LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - turbulent_adv_pdf.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 48 | 101 | 47.5 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 2 | 4 | 50.0 %

              Line data    Source code
           1             : !---------------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module turbulent_adv_pdf
           5             :
           6             :   ! Description:
           7             :   ! Calculates the turbulent advection term in the predictive equation for a
           8             :   ! variance or covariance where turbulent advection is calculated by
           9             :   ! integrating over the PDF.  This includes the following predictive fields:
          10             :   ! <w'thl'>, <w'rt'>, <rt'^2>, <thl'^2>, and <rt'thl'>, as well as passive
          11             :   ! scalar fields <w'sclr'>, <sclr'^2>, <sclr'rt'>, and <sclr'thl'>.  CLUBB
          12             :   ! does not produce a PDF for horizontal wind components u and v.  However, the
          13             :   ! horizontal wind variances <u'^2> and <v'^2> still use this code, as well.
          14             :
          15             :   ! References:
          16             :   !-------------------------------------------------------------------------
          17             :
          18             :   implicit none
          19             :
          20             :   public :: xpyp_term_ta_pdf_lhs,         &
          21             :             xpyp_term_ta_pdf_lhs_godunov, &
          22             :             xpyp_term_ta_pdf_rhs,         &
          23             :             xpyp_term_ta_pdf_rhs_godunov
          24             :
          25             :   private    ! Set default scope
          26             :
          27             :   integer, parameter :: &
          28             :     ndiags3 = 3
          29             :
          30             :   contains
          31             :
          32             :   !=============================================================================
          33    26805168 :   subroutine xpyp_term_ta_pdf_lhs( nz, ngrdcol, gr, coef_wpxpyp_implicit, & ! In
          34    26805168 :                                         rho_ds_zt, rho_ds_zm,            & ! In
          35    26805168 :                                         invrs_rho_ds_zm,                 & ! In
          36             :                                         l_upwind_xpyp_turbulent_adv,     & ! In
          37    26805168 :                                         sgn_turbulent_vel,               & ! In
          38    26805168 :                                         coef_wpxpyp_implicit_zm,         & ! In
          39    26805168 :                                         lhs_ta                           ) ! Out
          40             :
          41             :     ! Description:
          42             :     ! Turbulent advection of <w'x'>, <x'^2>, and <x'y'>:  implicit portion of
          43             :     ! the code.
          44             :     !
          45             :     ! 1) <w'x'>
          46             :     !
          47             :     ! The d<w'x'>/dt equation contains a turbulent advection term:
          48             :     !
          49             :     ! - (1/rho_ds) * d( rho_ds * <w'^2 x'> )/dz.
          50             :     !
          51             :     ! The value of <w'^2 x'> is found by integrating over the multivariate PDF
          52             :     ! of w and x, as detailed in function calc_wp2xp_pdf, which is found in
          53             :     ! module pdf_closure_module in pdf_closure_module.F90.
          54             :     !
          55             :     ! The equation obtained for <w'^2 x'> is written in terms of PDF parameters.
          56             :     ! Substitutions that are specific to the type of PDF used are made for the
          57             :     ! PDF parameters in order to write the <w'x'> turbulent advection term in
          58             :     ! the following form:
          59             :     !
          60             :     ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit.
          61             :     !
          62             :     ! For the ADG1 PDF, coef_wp2xp_implicit is a_1 * < w'^3 > / < w'^2 >, where
          63             :     ! a_1 is 1 / ( 1 - sigma_sqd_w ).  The value of term_wp2xp_explicit is 0, as
          64             :     ! the <w'x'> turbulent advection term is entirely implicit.
          65             :     !
          66             :     ! For the new PDF, the calculations of both coef_wp2xp_implicit and
          67             :     ! term_wp2xp_explicit are detailed in function calc_coefs_wp2xp_semiimpl,
          68             :     ! which is found in module new_pdf in new_pdf.F90.
          69             :     !
          70             :     ! For the new hybrid PDF, the calculation of coef_wp2xp_implicit is
          71             :     ! detailed in function calc_coef_wp2xp_implicit, which is found in module
          72             :     ! new_hybrid_pdf in new_hybrid_pdf.F90.  The value of term_wp2xp_explicit
          73             :     ! is 0, as the <w'x'> turbulent advection term is entirely implicit.
          74             :     !
          75             :     ! For explicit turbulent advection, the value of coef_wp2xp_implicit is 0
          76             :     ! and the value of term_wp2xp_explicit is <w'^2 x'>, as calculated by
          77             :     ! retaining the equation for <w'^2 x'> that is written in terms of PDF
          78             :     ! parameters.  This is a general form that can be used with any type of PDF.
          79             :     !
          80             :     ! The <w'x'> turbulent advection term is rewritten as:
          81             :     !
          82             :     ! - (1/rho_ds)
          83             :     !   * d( rho_ds * ( coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit ) )
          84             :     !     /dz.
          85             :     !
          86             :     ! The variable <w'x'> is evaluated at the (t+1) timestep, which allows the
          87             :     ! <w'x'> turbulent advection term to be expressed semi-implicitly as:
          88             :     !
          89             :     ! - (1/rho_ds) * d( rho_ds * coef_wp2xp_implicit * <w'x'>(t+1) )/dz
          90             :     ! - (1/rho_ds) * d( rho_ds * term_wp2xp_explicit )/dz.
          91             :     !
          92             :     ! The implicit portion of <w'x'> turbulent advection term is:
          93             :     !
          94             :     ! - (1/rho_ds) * d( rho_ds * coef_wp2xp_implicit * <w'x'>(t+1) )/dz.
          95             :     !
          96             :     ! Note:  When the implicit term is brought over to the left-hand side, the
          97             :     !        sign is reversed and the leading "-" in front of the implicit
          98             :     !        d[ ] / dz term is changed to a "+".
          99             :     !
         100             :     ! The timestep index (t+1) means that the value of <w'x'> being used is from
         101             :     ! the next timestep, which is being advanced to in solving the d<w'x'>/dt
         102             :     ! equation.
         103             :     !
         104             :     ! 2) <x'^2>
         105             :     !
         106             :     ! The d<x'^2>/dt equation contains a turbulent advection term:
         107             :     !
         108             :     ! - (1/rho_ds) * d( rho_ds * <w'x'^2> )/dz;
         109             :     !
         110             :     ! The value of <w'x'^2> is found by integrating over the multivariate PDF of
         111             :     ! w and x, as detailed in function calc_wpxp2_pdf, which is found in module
         112             :     ! pdf_closure_module in pdf_closure_module.F90.
         113             :     !
         114             :     ! The equation obtained for <w'x'^2> is written in terms of PDF parameters.
         115             :     ! Substitutions that are specific to the type of PDF used are made for the
         116             :     ! PDF parameters in order to write the <x'^2> turbulent advection term in
         117             :     ! the following form:
         118             :     !
         119             :     ! <w'x'^2> = coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit.
         120             :     !
         121             :     ! For the ADG1 PDF, the value of coef_wpxp2_implicit is
         122             :     ! (1/3)*beta * a_1 * < w'^3 > / < w'^2 >.  The value of term_wpxp2_explicit
         123             :     ! is (1-(1/3)*beta) * a_1^2 * < w'x' >^2 * < w'^3 > / < w'^2 >^2, where
         124             :     ! a_1 is 1 / ( 1 - sigma_sqd_w ).
         125             :     !
         126             :     ! For the new PDF, the calculation of coef_wpxp2_implicit is detailed in
         127             :     ! function calc_coef_wpxp2_implicit, which is found in module new_pdf in
         128             :     ! new_pdf.F90.  The value of term_wpxp2_explicit is 0, as the <x'^2>
         129             :     ! turbulent advection term is entirely implicit.
         130             :     !
         131             :     ! For the new hybrid PDF, the calculation of both coef_wpxp2_implicit and
         132             :     ! term_wpxp2_explicit are detailed in subroutine calc_coefs_wpxp2_semiimpl,
         133             :     ! which is found in module new_hybrid_pdf in new_hybrid_pdf.F90.
         134             :     !
         135             :     ! For explicit turbulent advection, the value of coef_wpxp2_implicit is 0
         136             :     ! and the value of term_wpxp2_explicit is <w'x'^2>, as calculated by
         137             :     ! retaining the equation for <w'x'^2> that is written in terms of PDF
         138             :     ! parameters.  This is a general form that can be used with any type of PDF.
         139             :     !
         140             :     ! The <x'^2> turbulent advection term is rewritten as:
         141             :     !
         142             :     ! - (1/rho_ds)
         143             :     !   * d( rho_ds * ( coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit ) )
         144             :     !     /dz;
         145             :     !
         146             :     ! The variable <x'^2> is evaluated at the (t+1) timestep, which allows the
         147             :     ! <x'^2> turbulent advection term to be expressed semi-implicitly as:
         148             :     !
         149             :     ! - (1/rho_ds) * d( rho_ds * coef_wpxp2_implicit * <x'^2>(t+1) )/dz;
         150             :     ! - (1/rho_ds) * d( rho_ds * term_wpxp2_explicit )/dz.
         151             :     !
         152             :     ! The implicit portion of <x'^2> turbulent advection term is:
         153             :     !
         154             :     ! - (1/rho_ds) * d( rho_ds * coef_wpxp2_implicit * <x'^2>(t+1) )/dz.
         155             :     !
         156             :     ! Note:  When the implicit term is brought over to the left-hand side, the
         157             :     !        sign is reversed and the leading "-" in front of all implicit
         158             :     !        d[ ] / dz terms is changed to a "+".
         159             :     !
         160             :     ! The timestep index (t+1) means that the value of <x'^2> being used is from
         161             :     ! the next timestep, which is being advanced to in solving the d<x'^2>/dt
         162             :     ! equation.
         163             :     !
         164             :     ! 3) <x'y'>
         165             :     !
         166             :     ! The d<x'y'>/dt equation contains a turbulent advection term:
         167             :     !
         168             :     ! - (1/rho_ds) * d( rho_ds * <w'x'y'> )/dz.
         169             :     !
         170             :     ! The value of <w'x'y'> is found by integrating over the multivariate PDF of
         171             :     ! w, x, and y, as detailed in function calc_wpxpyp_pdf, which is found in
         172             :     ! module pdf_closure_module in pdf_closure_module.F90.
         173             :     !
         174             :     ! The equation obtained for <w'x'y'> is written in terms of PDF parameters.
         175             :     ! Substitutions that are specific to the type of PDF used are made for the
         176             :     ! PDF parameters in order to write the <x'y'> turbulent advection term in
         177             :     ! the following form:
         178             :     !
         179             :     ! <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit.
         180             :     !
         181             :     ! For the ADG1 PDF, the value of coef_wpxpyp_implicit is
         182             :     ! (1/3)*beta * a_1 * < w'^3 > / < w'^2 >.  The value of term_wpxpyp_explicit
         183             :     ! is (1-(1/3)*beta) * a_1^2 * < w'x' > * < w'y' > * < w'^3 > / < w'^2 >^2,
         184             :     ! where a_1 is 1 / ( 1 - sigma_sqd_w ).
         185             :     !
         186             :     ! For the new PDF, the calculation of both coef_wpxpyp_implicit and
         187             :     ! term_wpxpyp_explicit are detailed in function calc_coefs_wpxpyp_semiimpl,
         188             :     ! which is found in module new_pdf in new_pdf.F90.
         189             :     !
         190             :     ! For the new hybrid PDF, the calculation of both coef_wpxpyp_implicit
         191             :     ! and term_wpxpyp_explicit are detailed in subroutine
         192             :     ! calc_coefs_wpxpyp_semiimpl, which is found in module new_hybrid_pdf in
         193             :     ! new_hybrid_pdf.F90.
         194             :     !
         195             :     ! For explicit turbulent advection, the value of coef_wpxpyp_implicit is 0
         196             :     ! and the value of term_wpxpyp_explicit is <w'x'y'>, as calculated by
         197             :     ! retaining the equation for <w'x'y'> that is written in terms of PDF
         198             :     ! parameters.  This is a general form that can be used with any type of PDF.
         199             :     !
         200             :     ! The <x'y'> turbulent advection term is rewritten as:
         201             :     !
         202             :     ! - (1/rho_ds)
         203             :     !   * d( rho_ds * ( coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit ) )
         204             :     !     /dz.
         205             :     !
         206             :     ! The variable <x'y'> is evaluated at the (t+1) timestep, which allows the
         207             :     ! <x'y'> turbulent advection term to be expressed semi-implicitly as:
         208             :     !
         209             :     ! - (1/rho_ds) * d( rho_ds * coef_wpxpyp_implicit * <x'y'>(t+1) )/dz
         210             :     ! - (1/rho_ds) * d( rho_ds * term_wpxpyp_explicit )/dz.
         211             :     !
         212             :     ! The implicit portion of <x'y'> turbulent advection term is:
         213             :     !
         214             :     ! - (1/rho_ds) * d( rho_ds * coef_wpxpyp_implicit * <x'y'>(t+1) )/dz.
         215             :     !
         216             :     ! Note:  When the implicit term is brought over to the left-hand side, the
         217             :     !        sign is reversed and the leading "-" in front of all implicit
         218             :     !        d[ ] / dz terms is changed to a "+".
         219             :     !
         220             :     ! The timestep index (t+1) means that the value of <x'y'> being used is from
         221             :     ! the next timestep, which is being advanced to in solving the d<x'y'>/dt
         222             :     ! equation.
         223             :     !
         224             :     ! When x and y are the same variable, <x'y'> reduces to <x'^2> and <w'x'y'>
         225             :     ! reduces to <w'x'^2>.  Likewise, when y is set equal to w, <x'y'> becomes
         226             :     ! <w'x'> and <w'x'y'> reduces to <w'^2 x'>.  The discretization and the code
         227             :     ! used in this function will be written generally in terms of <x'y'> and
         228             :     ! coef_wpxpyp_implicit, but also applies to <x'^2> and coef_wpxp2_implicit,
         229             :     ! as well as to <w'x'> and coef_wp2xp_implicit.
         230             :     !
         231             :     ! The implicit discretization of this term is as follows:
         232             :     !
         233             :     ! 1) Centered Discretization
         234             :     !
         235             :     ! The values of <x'y'> are found on the momentum levels, while the values of
         236             :     ! coef_wpxpyp_implicit are found on the thermodynamic levels, which is where
         237             :     ! they were originally calculated by the PDF.  Additionally, the values of
         238             :     ! rho_ds_zt are found on the thermodynamic levels, and the values of
         239             :     ! invrs_rho_ds_zm are found on the momentum levels.  The values of <x'y'>
         240             :     ! are interpolated to the intermediate thermodynamic levels as <x'y'>|_zt.
         241             :     ! At the thermodynamic levels, the values of coef_wpxpyp_implicit are
         242             :     ! multiplied by <x'y'>|_zt, and their product is multiplied by rho_ds_zt.
         243             :     ! Then, the derivative (d/dz) of that expression is taken over the central
         244             :     ! momentum level, where it is multiplied by -invrs_rho_ds_zm.  This yields
         245             :     ! the desired result.
         246             :     !
         247             :     ! =xpyp============================================================== m(k+1)
         248             :     !
         249             :     ! -xpyp_zt(interp)-------coef_wpxpyp_implicit---------rho_ds_zt------ t(k+1)
         250             :     !
         251             :     ! =xpyp=d(rho_ds_zt*coef_wpxpyp_implicit*xpyp_zt)/dz=invrs_rho_ds_zm= m(k)
         252             :     !
         253             :     ! -xpyp_zt(interp)-------coef_wpxpyp_implicit---------rho_ds_zt------ t(k)
         254             :     !
         255             :     ! =xpyp============================================================== m(k-1)
         256             :     !
         257             :     ! The vertical indices m(k+1), t(k+1), m(k), t(k), and m(k-1) correspond
         258             :     ! with altitudes zm(k+1), zt(k+1), zm(k), zt(k), and zm(k-1), respectively.
         259             :     ! The letter "t" is used for thermodynamic levels and the letter "m" is used
         260             :     ! for momentum levels.
         261             :     !
         262             :     ! invrs_dzm(k) = 1 / ( zt(k+1) - zt(k) )
         263             :     !
         264             :     ! 2) "Upwind" Discretization.
         265             :     !
         266             :     ! The values of <x'y'> are found on the momentum levels.  The values of
         267             :     ! coef_wpxpyp_implicit are originally calculated by the PDF on the
         268             :     ! thermodynamic levels.  They are interpolated to the intermediate momentum
         269             :     ! levels as coef_wpxpyp_implicit_zm.  Additionally, the values of rho_ds_zm
         270             :     ! and the values of invrs_rho_ds_zm are found on the momentum levels.  The
         271             :     ! sign of the turbulent velocity is found on the central momentum level.  At
         272             :     ! the momentum levels, the values of coef_wpxpyp_implicit_zm are multiplied
         273             :     ! by <x'y'>, and their product is multiplied by rho_ds_zm.  Then, the
         274             :     ! derivative (d/dz) of that expression is taken.  When the sign of the
         275             :     ! turbulent velocity is positive, the "wind" is coming from below, and the
         276             :     ! derivative involves the central momentum level and the momentum level
         277             :     ! immediately below it.  When the sign of the turbulent velocity is
         278             :     ! negative, the "wind" is coming from above, and the derivative involves the
         279             :     ! central momentum level and the momenum level immediately above it.  After
         280             :     ! the derivative is taken, it is multiplied by -invrs_rho_ds_zm at the
         281             :     ! central momentum level.  This yields the desired result.
         282             :     !
         283             :     ! The turbulent velocity for <x'y'> is <w'x'y'> / <x'y'>, which has units of
         284             :     ! m/s.  The sign of the turbulent velocity is sgn( <w'x'y'> / <x'y'> ),
         285             :     ! where:
         286             :     !
         287             :     ! sgn(x) = | 1; when x >= 0
         288             :     !          | -1; when x < 0.
         289             :     !
         290             :     ! The sign of the turbulent velocity can also be rewritten as
         291             :     ! sgn( <w'x'y'> ) / sgn( <x'y'> ).  When a variance (<x'^2>) is being solved
         292             :     ! for, y = x, and sgn( <x'^2> ) is always 1.  The sign of the turbulent
         293             :     ! velocity reduces to simply sgn( <w'x'^2> ).
         294             :     !
         295             :     ! ---------coef_wpxpyp_implicit-------------------------------------- t(k+2)
         296             :     !
         297             :     ! =xpyp=======coef_wpxpyp_implicit_zm(interp.)=======rho_ds_zm======= m(k+1)
         298             :     !
         299             :     ! ---------coef_wpxpyp_implicit-------------------------------------- t(k+1)
         300             :     !
         301             :     ! =xpyp===coef_wpxpyp_implicit_zm(interp.)=rho_ds_zm=invrs_rho_ds_zm= m(k)
         302             :     !
         303             :     ! ---------coef_wpxpyp_implicit-------------------------------------- t(k)
         304             :     !
         305             :     ! =xpyp=======coef_wpxpyp_implicit_zm(interp.)=======rho_ds_zm======= m(k-1)
         306             :     !
         307             :     ! ---------coef_wpxpyp_implicit-------------------------------------- t(k-1)
         308             :     !
         309             :     ! The vertical indices t(k+2), m(k+1), t(k+1), m(k), t(k), m(k-1), and
         310             :     ! t(k-1) correspond with altitudes zt(k+2), zm(k+1), zt(k+1), zm(k), zt(k),
         311             :     ! zm(k-1), and zt(k-1), respectively.  The letter "t" is used for
         312             :     ! thermodynamic levels and the letter "m" is used for momentum levels.
         313             :     !
         314             :     ! invrs_dzt(k+1) = 1 / ( zm(k+1) - zm(k) ); and
         315             :     ! invrs_dzt(k) = 1 / ( zm(k) - zm(k-1) ).
         316             :
         317             :     ! References:
         318             :     !-----------------------------------------------------------------------
         319             :
         320             :     use grid_class, only:  & ! gr%weights_zm2zt
         321             :         grid ! Type
         322             :
         323             :     use constants_clubb, only: &
         324             :         zero    ! Variable(s)
         325             :
         326             :     use clubb_precision, only: &
         327             :         core_rknd    ! Variable(s)
         328             :
         329             :     implicit none
         330             :
         331             :     ! Constant parameters
         332             :     integer, parameter :: &
         333             :       kp1_mdiag = 1, & ! Momentum superdiagonal index.
         334             :       k_mdiag   = 2, & ! Momentum main diagonal index.
         335             :       km1_mdiag = 3    ! Momentum subdiagonal index.
         336             :
         337             :     integer, parameter :: &
         338             :       m_above = 1, & ! Index for upper momentum level grid weight.
         339             :       m_below = 2    ! Index for lower momentum level grid weight.
         340             :
         341             :     ! ------------------------------ Input Variables ------------------------------
         342             :     integer, intent(in) :: &
         343             :       nz, &
         344             :       ngrdcol
         345             :
         346             :     type (grid), target, intent(in) :: gr
         347             :
         348             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         349             :       coef_wpxpyp_implicit,   & ! Coef. of <x'y'> in <w'x'y'>; t-levs  [m/s]
         350             :       rho_ds_zt,              & ! Dry, static density at t-levels      [kg/m^3]
         351             :       invrs_rho_ds_zm           ! Inv dry, static density at m-levels  [m^3/kg]
         352             :
         353             :     logical, intent(in) :: &
         354             :       l_upwind_xpyp_turbulent_adv    ! Flag to use "upwind" discretization
         355             :
         356             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         357             :       sgn_turbulent_vel,       & ! Sign of the turbulent velocity      [-]
         358             :       coef_wpxpyp_implicit_zm, & ! coef_wpxpyp_implicit interp m-levs  [m/s]
         359             :       rho_ds_zm                  ! Dry, static density at m-levs       [kg/m^3]
         360             :
         361             :     ! ------------------------------ Return Variable ------------------------------
         362             :     real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz), intent(out) :: &
         363             :       lhs_ta    ! LHS coefficient of xpyp turbulent advection  [1/s]
         364             :
         365             :     ! ------------------------------ Local Variables ------------------------------
         366             :     integer :: i, k, b    ! Vertical level index
         367             :
         368             :     ! ------------------------------ Begin Code ------------------------------
         369             :
         370             :     !$acc data copyin( gr, gr%weights_zm2zt, gr%invrs_dzm, gr%invrs_dzt, &
         371             :     !$acc              coef_wpxpyp_implicit, rho_ds_zt, invrs_rho_ds_zm, &
         372             :     !$acc              sgn_turbulent_vel, coef_wpxpyp_implicit_zm, rho_ds_zm ) &
         373             :     !$acc      copyout( lhs_ta )
         374             :
         375             :
         376             :     ! Set lower boundary array to 0
         377             :     !$acc parallel loop gang vector collapse(2) default(present)
         378   447583968 :     do i = 1, ngrdcol
         379  1709920368 :       do b = 1, ndiags3
         380  1683115200 :         lhs_ta(b,i,1) = zero
         381             :       end do
         382             :     end do
         383             :     !$acc end parallel loop
         384             :
         385    26805168 :     if ( .not. l_upwind_xpyp_turbulent_adv ) then
         386             :
         387             :       ! Centered discretization.
         388             :       !$acc parallel loop gang vector collapse(2) default(present)
         389   750544704 :       do k = 2, nz-1, 1
         390 12392091504 :         do i = 1, ngrdcol
         391             :
         392             :           ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
         393 23283093600 :           lhs_ta(kp1_mdiag,i,k) &
         394           0 :           = invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
         395 11641546800 :             * rho_ds_zt(i,k+1) * coef_wpxpyp_implicit(i,k+1) &
         396 46566187200 :             * gr%weights_zm2zt(i,k+1,m_above)
         397             :
         398             :           ! Momentum main diagonal: [ x xpyp(k,<t+1>) ]
         399             :           lhs_ta(k_mdiag,i,k) &
         400           0 :           = invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
         401             :             * ( rho_ds_zt(i,k+1) * coef_wpxpyp_implicit(i,k+1) &
         402           0 :                 * gr%weights_zm2zt(i,k+1,m_below) &
         403             :                 - rho_ds_zt(i,k) * coef_wpxpyp_implicit(i,k) &
         404 11641546800 :                   * gr%weights_zm2zt(i,k,m_above) )
         405             :
         406             :           ! Momentum subdiagonal: [ x xpyp(k-1,<t+1>) ]
         407             :           lhs_ta(km1_mdiag,i,k) &
         408           0 :           = - invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
         409             :               * rho_ds_zt(i,k) * coef_wpxpyp_implicit(i,k) &
         410 12383156448 :               * gr%weights_zm2zt(i,k,m_below)
         411             :
         412             :         end do
         413             :       end do
         414             :       !$acc end parallel loop
         415             :
         416             :     else ! l_upwind_xpyp_turbulent_adv
         417             :
         418             :       ! "Upwind" discretization
         419             :       !$acc parallel loop gang vector collapse(2) default(present)
         420  1501089408 :       do k = 2, nz-1, 1
         421 24784183008 :         do i = 1, ngrdcol
         422             :
         423 24766312896 :           if ( sgn_turbulent_vel(i,k) > zero ) then
         424             :
         425             :             ! The "wind" is blowing upward.
         426             :
         427             :             ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
         428  4946413726 :             lhs_ta(kp1_mdiag,i,k) = zero
         429             :
         430             :             ! Momentum main diagonal: [ x xpyp(k,<t+1>) ]
         431             :             lhs_ta(k_mdiag,i,k) &
         432           0 :              = invrs_rho_ds_zm(i,k) * gr%invrs_dzt(i,k) &
         433  4946413726 :                * rho_ds_zm(i,k) * coef_wpxpyp_implicit_zm(i,k)
         434             :
         435             :             ! Momentum subdiagonal: [ x xpyp(k-1,<t+1>) ]
         436             :             lhs_ta(km1_mdiag,i,k) &
         437           0 :              = - invrs_rho_ds_zm(i,k) * gr%invrs_dzt(i,k) &
         438  4946413726 :                  * rho_ds_zm(i,k-1) * coef_wpxpyp_implicit_zm(i,k-1)
         439             :
         440             :           else ! sgn_turbulent_vel < 0
         441             :
         442             :             ! The "wind" is blowing downward.
         443             :
         444             :             ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
         445             :             lhs_ta(kp1_mdiag,i,k) &
         446           0 :              = invrs_rho_ds_zm(i,k) * gr%invrs_dzt(i,k+1) &
         447 18336679874 :                * rho_ds_zm(i,k+1) * coef_wpxpyp_implicit_zm(i,k+1)
         448             :
         449             :             ! Momentum main diagonal: [ x xpyp(k,<t+1>) ]
         450             :             lhs_ta(k_mdiag,i,k) &
         451           0 :              = - invrs_rho_ds_zm(i,k) * gr%invrs_dzt(i,k+1) &
         452 18336679874 :                  * rho_ds_zm(i,k) * coef_wpxpyp_implicit_zm(i,k)
         453             :
         454             :             ! Momentum subdiagonal: [ x xpyp(k-1,<t+1>) ]
         455 18336679874 :             lhs_ta(km1_mdiag,i,k) = zero
         456             :
         457             :           end if ! sgn_turbulent_vel
         458             :
         459             :         end do
         460             :       end do
         461             :       !$acc end parallel loop
         462             :
         463             :     endif
         464             :
         465             :     ! Set upper boundary array to
         466             :     !$acc parallel loop gang vector collapse(2) default(present)
         467   447583968 :     do i = 1, ngrdcol
         468  1709920368 :       do b = 1, ndiags3
         469  1683115200 :         lhs_ta(b,i,nz) = zero
         470             :       end do
         471             :     end do
         472             :     !$acc end parallel loop
         473             :
         474             :     !$acc end data
         475             :
         476    26805168 :     return
         477             :
         478             :   end subroutine xpyp_term_ta_pdf_lhs
         479             :
         480             :   !=============================================================================================
         481           0 :   subroutine xpyp_term_ta_pdf_lhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
         482           0 :                                                 coef_wpxpyp_implicit, & ! Intent(in)
         483           0 :                                                 invrs_rho_ds_zm, rho_ds_zm,  & ! Intent(in)
         484           0 :                                                 lhs_ta )
         485             :   ! Intent(out)
         486             :   ! Description:
         487             :   !   This subroutine is a revised version of xpyp_term_ta_pdf_lhs_all. The
         488             :   !   revisions are maded to use the  Godunov-like upwind scheme for the
         489             :   !   vertical discretization of the turbulent advection term. This subroutine
         490             :   !   returns an array of 3 dimensional arrays, one for every grid level not
         491             :   !   including
         492             :   !   boundary values.
         493             :   !
         494             :   ! Optional Arguements:
         495             :   !   The optional arguements can be used to override the default indices.
         496             :   !   from_level - low index, default 2
         497             :   !   to level   - high index, default gr%nz-1
         498             :   !
         499             :   ! Notes:
         500             :   !   This subroutine exists for testing of Godunov-like upwind scheme.
         501             :   !   THIS SUBROUTINE DOES NOT HANDLE BOUNDARY CONDITIONS AND SETS THEM TO 0
         502             :   !---------------------------------------------------------------------------------------------
         503             :
         504             :     use grid_class, only:  & ! for gr%weights_zm2zt
         505             :     grid ! Type
         506             :
         507             :     use clubb_precision, only: &
         508             :         core_rknd    ! Variable(s)
         509             :
         510             :     implicit none
         511             :
         512             :     integer, parameter :: &
         513             :       kp1_mdiag = 1, & ! Momentum superdiagonal index.
         514             :       k_mdiag   = 2, & ! Momentum main diagonal index.
         515             :       km1_mdiag = 3    ! Momentum subdiagonal index.
         516             :
         517             :     !------------------- Input Variables -------------------
         518             :     integer, intent(in) :: &
         519             :       nz, &
         520             :       ngrdcol
         521             :
         522             :     type (grid), target, intent(in) :: gr
         523             :
         524             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         525             :         coef_wpxpyp_implicit,     & ! Coef. of <x'y'> in <w'x'y'>; t-lev [m/s]
         526             :         invrs_rho_ds_zm,          & ! Inv dry, static density @ m-level [m^3/kg]
         527             :         rho_ds_zm                   ! Dry, static density at m-lev [kg/m^3]
         528             :
         529             :     !------------------- Output Variables -------------------
         530             :     real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz), intent(out) :: &
         531             :         lhs_ta
         532             :
         533             :     !---------------- Local Variables -------------------
         534             :     integer :: &
         535             :         i, k, b          ! Loop variable for current grid level
         536             :
         537             :     !---------------- Begin Code -------------------
         538             :
         539             :     !$acc data copyin( gr, gr%invrs_dzm, &
         540             :     !$acc              coef_wpxpyp_implicit, invrs_rho_ds_zm, rho_ds_zm ) &
         541             :     !$acc      copyout( lhs_ta )
         542             :
         543             :     ! Set lower boundary array to 0
         544             :     !$acc parallel loop gang vector collapse(2) default(present)
         545           0 :     do i = 1, ngrdcol
         546           0 :       do b = 1, ndiags3
         547           0 :         lhs_ta(b,i,1) = 0.0_core_rknd
         548             :       end do
         549             :     end do
         550             :     !$acc end parallel loop
         551             :
         552             :     ! Godunov-like upwind discretization
         553             :     !$acc parallel loop gang vector collapse(2) default(present)
         554           0 :     do k = 2, nz-1
         555           0 :       do i = 1, ngrdcol
         556             :
         557             :         ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
         558           0 :         lhs_ta(kp1_mdiag,i,k) = invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
         559           0 :                                 * rho_ds_zm(i,k+1) &
         560           0 :                                 * min(0.0_core_rknd,coef_wpxpyp_implicit(i,k+1))
         561             :
         562             :         ! Momentum main diagonal: [ x xpyp(k,<t+1>) ]
         563           0 :         lhs_ta(k_mdiag,i,k) = invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
         564             :                               * rho_ds_zm(i,k) &
         565             :                               * ( max(0.0_core_rknd,coef_wpxpyp_implicit(i,k+1)) - &
         566           0 :                                   min(0.0_core_rknd,coef_wpxpyp_implicit(i,k)) )
         567             :
         568             :         ! Momentum subdiagonal: [ x xpyp(k-1,<t+1>) ]
         569           0 :         lhs_ta(km1_mdiag,i,k) = - invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
         570           0 :                                   * rho_ds_zm(i,k-1) &
         571           0 :                                   * max(0.0_core_rknd,coef_wpxpyp_implicit(i,k) )
         572             :       end do
         573             :     end do
         574             :     !$acc end parallel loop
         575             :
         576             :     ! Set upper boundary array to 0
         577             :     !$acc parallel loop gang vector collapse(2) default(present)
         578           0 :     do i = 1, ngrdcol
         579           0 :       do b = 1, ndiags3
         580           0 :         lhs_ta(b,i,nz) = 0.0_core_rknd
         581             :       end do
         582             :     end do
         583             :     !$acc end parallel loop
         584             :
         585             :     !$acc end data
         586             :
         587           0 :     return
         588             :
         589             :   end subroutine xpyp_term_ta_pdf_lhs_godunov
         590             :
         591             :   !=============================================================================
         592    44675280 :   subroutine xpyp_term_ta_pdf_rhs( nz, ngrdcol, gr, term_wpxpyp_explicit,  & ! In
         593    44675280 :                                         rho_ds_zt, rho_ds_zm,                   & ! In
         594    44675280 :                                         invrs_rho_ds_zm,                        & ! In
         595             :                                         l_upwind_xpyp_turbulent_adv,            & ! In
         596    44675280 :                                         sgn_turbulent_vel,                      & ! In
         597    44675280 :                                         term_wpxpyp_explicit_zm,                & ! In
         598    44675280 :                                         rhs_ta )                                  ! Out
         599             :
         600             :     ! Description:
         601             :     ! Turbulent advection of <w'x'>, <x'^2>, and <x'y'>:  explicit portion of
         602             :     ! the code.
         603             :     !
         604             :     ! 1) <w'x'>
         605             :     !
         606             :     ! The d<w'x'>/dt equation contains a turbulent advection term:
         607             :     !
         608             :     ! - (1/rho_ds) * d( rho_ds * <w'^2 x'> )/dz.
         609             :     !
         610             :     ! The value of <w'^2 x'> is found by integrating over the multivariate PDF
         611             :     ! of w and x, as detailed in function calc_wp2xp_pdf, which is found in
         612             :     ! module pdf_closure_module in pdf_closure_module.F90.
         613             :     !
         614             :     ! The equation obtained for <w'^2 x'> is written in terms of PDF parameters.
         615             :     ! Substitutions that are specific to the type of PDF used are made for the
         616             :     ! PDF parameters in order to write the <w'x'> turbulent advection term in
         617             :     ! the following form:
         618             :     !
         619             :     ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit.
         620             :     !
         621             :     ! For the ADG1 PDF, coef_wp2xp_implicit is a_1 * < w'^3 > / < w'^2 >, where
         622             :     ! a_1 is 1 / ( 1 - sigma_sqd_w ).  The value of term_wp2xp_explicit is 0, as
         623             :     ! the <w'x'> turbulent advection term is entirely implicit.
         624             :     !
         625             :     ! For the new PDF, the calculations of both coef_wp2xp_implicit and
         626             :     ! term_wp2xp_explicit are detailed in function calc_coefs_wp2xp_semiimpl,
         627             :     ! which is found in module new_pdf in new_pdf.F90.
         628             :     !
         629             :     ! For the new hybrid PDF, the calculation of coef_wp2xp_implicit is
         630             :     ! detailed in function calc_coef_wp2xp_implicit, which is found in module
         631             :     ! new_hybrid_pdf in new_hybrid_pdf.F90.  The value of term_wp2xp_explicit
         632             :     ! is 0, as the <w'x'> turbulent advection term is entirely implicit.
         633             :     !
         634             :     ! For explicit turbulent advection, the value of coef_wp2xp_implicit is 0
         635             :     ! and the value of term_wp2xp_explicit is <w'^2 x'>, as calculated by
         636             :     ! retaining the equation for <w'^2 x'> that is written in terms of PDF
         637             :     ! parameters.  This is a general form that can be used with any type of PDF.
         638             :     !
         639             :     ! The <w'x'> turbulent advection term is rewritten as:
         640             :     !
         641             :     ! - (1/rho_ds)
         642             :     !   * d( rho_ds * ( coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit ) )
         643             :     !     /dz.
         644             :     !
         645             :     ! The variable <w'x'> is evaluated at the (t+1) timestep, which allows the
         646             :     ! <w'x'> turbulent advection term to be expressed semi-implicitly as:
         647             :     !
         648             :     ! - (1/rho_ds) * d( rho_ds * coef_wp2xp_implicit * <w'x'>(t+1) )/dz
         649             :     ! - (1/rho_ds) * d( rho_ds * term_wp2xp_explicit )/dz.
         650             :     !
         651             :     ! The explicit portion of <w'x'> turbulent advection term is:
         652             :     !
         653             :     ! - (1/rho_ds) * d( rho_ds * term_wp2xp_explicit )/dz.
         654             :     !
         655             :     ! 2) <x'^2>
         656             :     !
         657             :     ! The d<x'^2>/dt equation contains a turbulent advection term:
         658             :     !
         659             :     ! - (1/rho_ds) * d( rho_ds * <w'x'^2> )/dz;
         660             :     !
         661             :     ! The value of <w'x'^2> is found by integrating over the multivariate PDF of
         662             :     ! w and x, as detailed in function calc_wpxp2_pdf, which is found in module
         663             :     ! pdf_closure_module in pdf_closure_module.F90.
         664             :     !
         665             :     ! The equation obtained for <w'x'^2> is written in terms of PDF parameters.
         666             :     ! Substitutions that are specific to the type of PDF used are made for the
         667             :     ! PDF parameters in order to write the <x'^2> turbulent advection term in
         668             :     ! the following form:
         669             :     !
         670             :     ! <w'x'^2> = coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit.
         671             :     !
         672             :     ! For the ADG1 PDF, the value of coef_wpxp2_implicit is
         673             :     ! (1/3)*beta * a_1 * < w'^3 > / < w'^2 >.  The value of term_wpxp2_explicit
         674             :     ! is (1-(1/3)*beta) * a_1^2 * < w'x' >^2 * < w'^3 > / < w'^2 >^2, where
         675             :     ! a_1 is 1 / ( 1 - sigma_sqd_w ).
         676             :     !
         677             :     ! For the new PDF, the calculation of coef_wpxp2_implicit is detailed in
         678             :     ! function calc_coef_wpxp2_implicit, which is found in module new_pdf in
         679             :     ! new_pdf.F90.  The value of term_wpxp2_explicit is 0, as the <x'^2>
         680             :     ! turbulent advection term is entirely implicit.
         681             :     !
         682             :     ! For the new hybrid PDF, the calculation of both coef_wpxp2_implicit and
         683             :     ! term_wpxp2_explicit are detailed in subroutine calc_coefs_wpxp2_semiimpl,
         684             :     ! which is found in module new_hybrid_pdf in new_hybrid_pdf.F90.
         685             :     !
         686             :     ! For explicit turbulent advection, the value of coef_wpxp2_implicit is 0
         687             :     ! and the value of term_wpxp2_explicit is <w'x'^2>, as calculated by
         688             :     ! retaining the equation for <w'x'^2> that is written in terms of PDF
         689             :     ! parameters.  This is a general form that can be used with any type of PDF.
         690             :     !
         691             :     ! The <x'^2> turbulent advection term is rewritten as:
         692             :     !
         693             :     ! - (1/rho_ds)
         694             :     !   * d( rho_ds * ( coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit ) )
         695             :     !     /dz;
         696             :     !
         697             :     ! The variable <x'^2> is evaluated at the (t+1) timestep, which allows the
         698             :     ! <x'^2> turbulent advection term to be expressed semi-implicitly as:
         699             :     !
         700             :     ! - (1/rho_ds) * d( rho_ds * coef_wpxp2_implicit * <x'^2>(t+1) )/dz;
         701             :     ! - (1/rho_ds) * d( rho_ds * term_wpxp2_explicit )/dz.
         702             :     !
         703             :     ! The explicit portion of <x'^2> turbulent advection term is:
         704             :     !
         705             :     ! - (1/rho_ds) * d( rho_ds * term_wpxp2_explicit )/dz.
         706             :     !
         707             :     ! 3) <x'y'>
         708             :     !
         709             :     ! The d<x'y'>/dt equation contains a turbulent advection term:
         710             :     !
         711             :     ! - (1/rho_ds) * d( rho_ds * <w'x'y'> )/dz.
         712             :     !
         713             :     ! The value of <w'x'y'> is found by integrating over the multivariate PDF of
         714             :     ! w, x, and y, as detailed in function calc_wpxpyp_pdf, which is found in
         715             :     ! module pdf_closure_module in pdf_closure_module.F90.
         716             :     !
         717             :     ! The equation obtained for <w'x'y'> is written in terms of PDF parameters.
         718             :     ! Substitutions that are specific to the type of PDF used are made for the
         719             :     ! PDF parameters in order to write the <x'y'> turbulent advection term in
         720             :     ! the following form:
         721             :     !
         722             :     ! <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit.
         723             :     !
         724             :     ! For the ADG1 PDF, the value of coef_wpxpyp_implicit is
         725             :     ! (1/3)*beta * a_1 * < w'^3 > / < w'^2 >.  The value of term_wpxpyp_explicit
         726             :     ! is (1-(1/3)*beta) * a_1^2 * < w'x' > * < w'y' > * < w'^3 > / < w'^2 >^2,
         727             :     ! where a_1 is 1 / ( 1 - sigma_sqd_w ).
         728             :     !
         729             :     ! For the new PDF, the calculation of both coef_wpxpyp_implicit and
         730             :     ! term_wpxpyp_explicit are detailed in function calc_coefs_wpxpyp_semiimpl,
         731             :     ! which is found in module new_pdf in new_pdf.F90.
         732             :     !
         733             :     ! For the new hybrid PDF, the calculation of both coef_wpxpyp_implicit
         734             :     ! and term_wpxpyp_explicit are detailed in subroutine
         735             :     ! calc_coefs_wpxpyp_semiimpl, which is found in module new_hybrid_pdf in
         736             :     ! new_hybrid_pdf.F90.
         737             :     !
         738             :     ! For explicit turbulent advection, the value of coef_wpxpyp_implicit is 0
         739             :     ! and the value of term_wpxpyp_explicit is <w'x'y'>, as calculated by
         740             :     ! retaining the equation for <w'x'y'> that is written in terms of PDF
         741             :     ! parameters.  This is a general form that can be used with any type of PDF.
         742             :     !
         743             :     ! The <x'y'> turbulent advection term is rewritten as:
         744             :     !
         745             :     ! - (1/rho_ds)
         746             :     !   * d( rho_ds * ( coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit ) )
         747             :     !     /dz.
         748             :     !
         749             :     ! The variable <x'y'> is evaluated at the (t+1) timestep, which allows the
         750             :     ! <x'y'> turbulent advection term to be expressed semi-implicitly as:
         751             :     !
         752             :     ! - (1/rho_ds) * d( rho_ds * coef_wpxpyp_implicit * <x'y'>(t+1) )/dz
         753             :     ! - (1/rho_ds) * d( rho_ds * term_wpxpyp_explicit )/dz.
         754             :     !
         755             :     ! The explicit portion of <x'y'> turbulent advection term is:
         756             :     !
         757             :     ! - (1/rho_ds) * d( rho_ds * term_wpxpyp_explicit )/dz.
         758             :     !
         759             :     ! When x and y are the same variable, <x'y'> reduces to <x'^2> and <w'x'y'>
         760             :     ! reduces to <w'x'^2>.  Likewise, when y is set equal to w, <x'y'> becomes
         761             :     ! <w'x'> and <w'x'y'> reduces to <w'^2 x'>.  The discretization and the code
         762             :     ! used in this function will be written generally in terms of <x'y'> and
         763             :     ! term_wpxpyp_explicit, but also applies to <x'^2> and term_wpxp2_explicit,
         764             :     ! as well as to <w'x'> and term_wp2xp_explicit.
         765             :     !
         766             :     ! The explicit discretization of this term is as follows:
         767             :     !
         768             :     ! 1) Centered Discretization
         769             :     !
         770             :     ! The values of <x'y'> are found on the momentum levels, while the values of
         771             :     ! term_wpxpyp_explicit are found on the thermodynamic levels, which is where
         772             :     ! they were originally calculated by the PDF.  Additionally, the values of
         773             :     ! rho_ds_zt are found on the thermodynamic levels, and the values of
         774             :     ! invrs_rho_ds_zm are found on the momentum levels.  At the thermodynamic
         775             :     ! levels, the values of term_wpxpyp_explicit are multiplied by rho_ds_zt.
         776             :     ! Then, the derivative (d/dz) of that expression is taken over the central
         777             :     ! momentum level, where it is multiplied by -invrs_rho_ds_zm.  This yields
         778             :     ! the desired result.
         779             :     !
         780             :     ! ---------term_wpxpyp_explicitp1-------rho_ds_ztp1------------------ t(k+1)
         781             :     !
         782             :     ! ==xpyp==d( rho_ds_zt * term_wpxpyp_explicit )/dz==invrs_rho_ds_zm== m(k)
         783             :     !
         784             :     ! ---------term_wpxpyp_explicit---------rho_ds_zt-------------------- t(k)
         785             :     !
         786             :     ! The vertical indices t(k+1), m(k), and t(k) correspond with altitudes
         787             :     ! zt(k+1), zm(k), and zt(k), respectively.  The letter "t" is used for
         788             :     ! thermodynamic levels and the letter "m" is used for momentum levels.
         789             :     !
         790             :     ! invrs_dzm(k) = 1 / ( zt(k+1) - zt(k) )
         791             :     !
         792             :     ! 2) "Upwind" Discretization.
         793             :     !
         794             :     ! The values of <x'y'> are found on the momentum levels.  The values of
         795             :     ! term_wpxpyp_explicit are originally calculated by the PDF on the
         796             :     ! thermodynamic levels.  They are interpolated to the intermediate momentum
         797             :     ! levels as term_wpxpyp_explicit_zm.  Additionally, the values of rho_ds_zm
         798             :     ! and the values of invrs_rho_ds_zm are found on the momentum levels.  The
         799             :     ! sign of the turbulent velocity is found on the central momentum level.  At
         800             :     ! the momentum levels, the values of term_wpxpyp_explicit_zm are multiplied
         801             :     ! by rho_ds_zm.  Then, the derivative (d/dz) of that expression is taken.
         802             :     ! When the sign of the turbulent velocity is positive, the "wind" is coming
         803             :     ! from below, and the derivative involves the central momentum level and the
         804             :     ! momentum level immediately below it.  When the sign of the turbulent
         805             :     ! velocity is negative, the "wind" is coming from above, and the derivative
         806             :     ! involves the central momentum level and the momenum level immediately
         807             :     ! above it.  After the derivative is taken, it is multiplied by
         808             :     ! -invrs_rho_ds_zm at the central momentum level.  This yields the desired
         809             :     ! result.
         810             :     !
         811             :     ! The turbulent velocity for <x'y'> is <w'x'y'> / <x'y'>, which has units of
         812             :     ! m/s.  The sign of the turbulent velocity is sgn( <w'x'y'> / <x'y'> ),
         813             :     ! where:
         814             :     !
         815             :     ! sgn(x) = | 1; when x >= 0
         816             :     !          | -1; when x < 0.
         817             :     !
         818             :     ! The sign of the turbulent velocity can also be rewritten as
         819             :     ! sgn( <w'x'y'> ) / sgn( <x'y'> ).  When a variance (<x'^2>) is being solved
         820             :     ! for, y = x, and sgn( <x'^2> ) is always 1.  The sign of the turbulent
         821             :     ! velocity reduces to simply sgn( <w'x'^2> ).
         822             :     !
         823             :     ! ---------term_wpxpyp_explicit-------------------------------------- t(k+2)
         824             :     !
         825             :     ! ========term_wpxpyp_explicit_zm(interp.)=======rho_ds_zm=========== m(k+1)
         826             :     !
         827             :     ! ---------term_wpxpyp_explicit-------------------------------------- t(k+1)
         828             :     !
         829             :     ! =xpyp===term_wpxpyp_explicit_zm(interp.)=rho_ds_zm=invrs_rho_ds_zm= m(k)
         830             :     !
         831             :     ! ---------term_wpxpyp_explicit-------------------------------------- t(k)
         832             :     !
         833             :     ! ========term_wpxpyp_explicit_zm(interp.)=======rho_ds_zm=========== m(k-1)
         834             :     !
         835             :     ! ---------term_wpxpyp_explicit-------------------------------------- t(k-1)
         836             :     !
         837             :     ! The vertical indices t(k+2), m(k+1), t(k+1), m(k), t(k), m(k-1), and
         838             :     ! t(k-1) correspond with altitudes zt(k+2), zm(k+1), zt(k+1), zm(k), zt(k),
         839             :     ! zm(k-1), and zt(k-1), respectively.  The letter "t" is used for
         840             :     ! thermodynamic levels and the letter "m" is used for momentum levels.
         841             :     !
         842             :     ! invrs_dzt(k+1) = 1 / ( zm(k+1) - zm(k) ); and
         843             :     ! invrs_dzt(k) = 1 / ( zm(k) - zm(k-1) ).
         844             :
         845             :     ! References:
         846             :     !-----------------------------------------------------------------------
         847             :
         848             :     use grid_class, only:  &
         849             :         grid ! Type
         850             :
         851             :     use constants_clubb, only: &
         852             :         zero    ! Variable(s)
         853             :
         854             :     use clubb_precision, only: &
         855             :         core_rknd ! Variable(s)
         856             :
         857             :     implicit none
         858             :
         859             :     ! ----------------------------- Input Variables -----------------------------
         860             :     integer, intent(in) :: &
         861             :       nz, &
         862             :       ngrdcol
         863             :
         864             :     type (grid), target, intent(in) :: gr
         865             :
         866             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         867             :       term_wpxpyp_explicit, & ! RHS: <w'x'y'> eq; t-levs      [m/s(x un)(y un)]
         868             :       rho_ds_zt,            & ! Dry, static density at t-levs          [kg/m^3]
         869             :       invrs_rho_ds_zm         ! Inv dry, static density at m-levs      [m^3/kg]
         870             :
         871             :     logical, intent(in) :: &
         872             :       l_upwind_xpyp_turbulent_adv    ! Flag to use "upwind" discretization
         873             :
         874             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         875             :       sgn_turbulent_vel,       & ! Sign of the turbulent velocity           [-]
         876             :       term_wpxpyp_explicit_zm, & ! term_wpxpyp_expl. zm       [m/s(x un)(y un)]
         877             :       rho_ds_zm                  ! Dry, static density at m-levs       [kg/m^3]
         878             :
         879             :     ! ----------------------------- Return Variable -----------------------------
         880             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         881             :       rhs_ta    ! RHS portion of xpyp turbulent advection      [(x un)(y un)/s]
         882             :
         883             :     ! ----------------------------- Local Variables -----------------------------
         884             :     integer :: i, k    ! Vertical level index
         885             :
         886             :     ! ----------------------------- Begin Code -----------------------------
         887             :
         888             :     !$acc data copyin( gr, gr%invrs_dzm, gr%invrs_dzt, &
         889             :     !$acc              term_wpxpyp_explicit, rho_ds_zt, invrs_rho_ds_zm, &
         890             :     !$acc              sgn_turbulent_vel, term_wpxpyp_explicit_zm, rho_ds_zm ) &
         891             :     !$acc      copyout( rhs_ta )
         892             :
         893             :
         894             :     ! Set lower boundary value to 0
         895             :     !$acc parallel loop gang vector default(present)
         896   745973280 :     do i = 1, ngrdcol
         897   745973280 :       rhs_ta(i,1) = zero
         898             :     end do
         899             :     !$acc end parallel loop
         900             :
         901    44675280 :     if ( .not. l_upwind_xpyp_turbulent_adv ) then
         902             :
         903             :       ! Centered discretization.
         904             :       !$acc parallel loop gang vector collapse(2) default(present)
         905           0 :       do k = 2, nz-1, 1
         906           0 :         do i = 1, ngrdcol
         907             :
         908           0 :           rhs_ta(i,k) &
         909             :           = - invrs_rho_ds_zm(i,k) &
         910           0 :               * gr%invrs_dzm(i,k) * ( rho_ds_zt(i,k+1) * term_wpxpyp_explicit(i,k+1) &
         911           0 :                                  - rho_ds_zt(i,k) * term_wpxpyp_explicit(i,k) )
         912             :         end do
         913             :       end do ! k = 2, nz-1, 1
         914             :       !$acc end parallel loop
         915             :
         916             :     else ! l_upwind_xpyp_turbulent_adv
         917             :
         918             :       ! "Upwind" discretization
         919             :       !$acc parallel loop gang vector collapse(2) default(present)
         920  3752723520 :       do k = 2, nz-1, 1
         921 61960457520 :         do i = 1, ngrdcol
         922             :
         923 61915782240 :           if ( sgn_turbulent_vel(i,k) > zero ) then
         924             :
         925             :              ! The "wind" is blowing upward.
         926             :              rhs_ta(i,k) &
         927             :              = - invrs_rho_ds_zm(i,k) &
         928           0 :                  * gr%invrs_dzt(i,k) &
         929             :                  * ( rho_ds_zm(i,k) * term_wpxpyp_explicit_zm(i,k) &
         930 12366034315 :                      - rho_ds_zm(i,k-1) * term_wpxpyp_explicit_zm(i,k-1) )
         931             :
         932             :           else ! sgn_turbulent_vel < 0
         933             :
         934             :              ! The "wind" is blowing downward.
         935             :              rhs_ta(i,k) &
         936             :              = - invrs_rho_ds_zm(i,k) &
         937           0 :                  * gr%invrs_dzt(i,k+1) &
         938 45841699685 :                  * ( rho_ds_zm(i,k+1) * term_wpxpyp_explicit_zm(i,k+1) &
         939 91683399370 :                      - rho_ds_zm(i,k) * term_wpxpyp_explicit_zm(i,k) )
         940             :
         941             :           endif ! sgn_turbulent_vel
         942             :
         943             :         end do
         944             :       end do ! k = 2, nz-1, 1
         945             :       !$acc end parallel loop
         946             :
         947             :     end if
         948             :
         949             :     ! Set upper boundary value to 0
         950             :     !$acc parallel loop gang vector default(present)
         951   745973280 :     do i = 1, ngrdcol
         952   745973280 :       rhs_ta(i,nz) = zero
         953             :     end do
         954             :     !$acc end parallel loop
         955             :
         956             :     !$acc end data
         957             :
         958    44675280 :     return
         959             :
         960             :   end subroutine xpyp_term_ta_pdf_rhs
         961             :
         962             :   !=============================================================================
         963           0 :   subroutine xpyp_term_ta_pdf_rhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
         964           0 :                                                 term_wpxpyp_explicit_zm, & ! Intent(in)
         965           0 :                                                 invrs_rho_ds_zm, & ! Intent(in)
         966           0 :                                                 sgn_turbulent_vel, & ! Intent(in)
         967           0 :                                                 rho_ds_zm, & ! Intent(in)
         968           0 :                                                 rhs_ta ) ! Intent(out)
         969             :     ! Description:
         970             :     !   This subroutine intends to add godunov upwind difference scheme based
         971             :     !   on xpyp_term_ta_pdf_rhs.  The revisions are maded to use the Godunov-like
         972             :     !   upwind scheme for the vertical discretization.
         973             :     !   This subroutine returns an array of values for every grid level.
         974             :     !
         975             :     ! Optional Arguements:
         976             :     !   The optional arguements can be used to override the default indices.
         977             :     !   from_level - low index, default 2
         978             :     !   to level - high index, default gr%nz-1
         979             :     !
         980             :     ! Notes:
         981             :     !   This subroutine exists for testing of Godunov-like upwind scheme.
         982             :     !   THIS SUBROUTINE DOES NOT HANDLE BOUNDARY CONDITIONS AND SETS THEM TO 0
         983             :     !--------------------------------------------------------------------------------------------
         984             :     use clubb_precision, only: &
         985             :       core_rknd ! Variable(s)
         986             :
         987             :     use grid_class, only: &
         988             :       grid ! Type
         989             :
         990             :     implicit none
         991             :
         992             :     !------------------- Input Variables -------------------
         993             :     integer, intent(in) :: &
         994             :       nz, &
         995             :       ngrdcol
         996             :
         997             :     type (grid), target, intent(in) :: gr
         998             :
         999             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1000             :       term_wpxpyp_explicit_zm,      & ! RHS: <w'x'y'> eq; m-lev(k)   [m/s(x un)(y un)]
        1001             :       invrs_rho_ds_zm,              & ! Inv dry, static density at m-lev (k) [m^3/kg]
        1002             :       sgn_turbulent_vel,            & ! Sign of the turbulent velocity [-]
        1003             :       rho_ds_zm                       ! Dry, static density at m-lev (k) [kg/m^3]
        1004             :
        1005             :     !------------------- Output Variables -------------------
        1006             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1007             :       rhs_ta
        1008             :
        1009             :     !---------------- Local Variables -------------------
        1010             :     integer :: &
        1011             :       i, k             ! Loop variable for current grid level
        1012             :
        1013             :     !---------------- Begin Code -------------------
        1014             :
        1015             :     !$acc data copyin( gr, gr%invrs_dzm, &
        1016             :     !$acc              term_wpxpyp_explicit_zm, invrs_rho_ds_zm, &
        1017             :     !$acc              sgn_turbulent_vel, rho_ds_zm ) &
        1018             :     !$acc      copyout( rhs_ta )
        1019             :
        1020             :     ! Set lower boundary value to 0
        1021             :     !$acc parallel loop gang vector default(present)
        1022           0 :     do i = 1, ngrdcol
        1023           0 :       rhs_ta(i,1) = 0.0_core_rknd
        1024             :     end do
        1025             :     !$acc end parallel loop
        1026             :
        1027             :     !$acc parallel loop gang vector collapse(2) default(present)
        1028           0 :     do k = 2, nz-1
        1029           0 :       do i = 1, ngrdcol
        1030           0 :         rhs_ta(i,k) = - invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
        1031           0 :                     * ( min(0.0_core_rknd, sgn_turbulent_vel(i,k+1)) &
        1032             :                           * rho_ds_zm(i,k+1) * term_wpxpyp_explicit_zm(i,k+1) &
        1033             :                       + max(0.0_core_rknd, sgn_turbulent_vel(i,k+1)) &
        1034             :                           * rho_ds_zm(i,k)   * term_wpxpyp_explicit_zm(i,k) &
        1035             :                       - min(0.0_core_rknd, sgn_turbulent_vel(i,k))&
        1036             :                           * rho_ds_zm(i,k)   * term_wpxpyp_explicit_zm(i,k) &
        1037             :                       - max(0.0_core_rknd, sgn_turbulent_vel(i,k)) &
        1038           0 :                           * rho_ds_zm(i,k-1) * term_wpxpyp_explicit_zm(i,k-1) &
        1039           0 :                       )
        1040             :       end do
        1041             :     end do
        1042             :     !$acc end parallel loop
        1043             :
        1044             :     ! Set upper boundary value to 0
        1045             :     !$acc parallel loop gang vector default(present)
        1046           0 :     do i = 1, ngrdcol
        1047           0 :       rhs_ta(i,nz) = 0.0_core_rknd
        1048             :     end do
        1049             :     !$acc end parallel loop
        1050             :
        1051             :     !$acc end data
        1052             :
        1053           0 :     return
        1054             :
        1055             :   end subroutine xpyp_term_ta_pdf_rhs_godunov
        1056             :
        1057             : !===============================================================================
        1058             :
        1059             : end module turbulent_adv_pdf

Generated by: LCOV version 1.14