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