LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - mean_adv.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 45 | 62 | 72.6 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 2 | 100.0 %

              Line data    Source code
           1             : !-----------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module mean_adv
           5             :
           6             :   ! Description:
           7             :   ! Module mean_adv computes the mean advection terms for all of the
           8             :   ! time-tendency (prognostic) equations in the CLUBB parameterization.  All of
           9             :   ! the mean advection terms are solved for completely implicitly, and therefore
          10             :   ! become part of the left-hand side of their respective equations.
          11             :   !
          12             :   ! Function term_ma_zt_lhs handles the mean advection terms for the variables
          13             :   ! located at thermodynamic grid levels.  These variables are:  rtm, thlm, wp3,
          14             :   ! all hydrometeor species, and sclrm.
          15             :   !
          16             :   ! Function term_ma_zm_lhs handles the mean advection terms for the variables
          17             :   ! located at momentum grid levels.  The variables are:  wprtp, wpthlp, wp2,
          18             :   ! rtp2, thlp2, rtpthlp, up2, vp2, wpsclrp, sclrprtp, sclrpthlp, and sclrp2.
          19             :
          20             :   implicit none
          21             :
          22             :   private ! Default scope
          23             :
          24             :   public :: term_ma_zt_lhs, &
          25             :             term_ma_zm_lhs
          26             :
          27             :   integer, parameter :: &
          28             :     ndiags3 = 3
          29             :
          30             :   contains
          31             :
          32             :   !=============================================================================
          33     8935056 :   subroutine term_ma_zt_lhs( nz, ngrdcol, wm_zt, weights_zt2zm, & ! Intent(in)
          34     8935056 :                                   invrs_dzt, invrs_dzm,     & ! Intent(in)
          35             :                                   l_upwind_xm_ma,           & ! Intent(in)
          36     8935056 :                                   lhs_ma )                    ! Intent(out)
          37             :
          38             :     ! Description:
          39             :     ! Mean advection of var_zt:  implicit portion of the code.
          40             :     !
          41             :     ! The variable "var_zt" stands for a variable that is located at
          42             :     ! thermodynamic grid levels.
          43             :     !
          44             :     ! The d(var_zt)/dt equation contains a mean advection term:
          45             :     !
          46             :     ! - w * d(var_zt)/dz.
          47             :     !
          48             :     ! This term is solved for completely implicitly, such that:
          49             :     !
          50             :     ! - w * d( var_zt(t+1) )/dz.
          51             :     !
          52             :     ! Note:  When the term is brought over to the left-hand side, the sign
          53             :     !        is reversed and the leading "-" in front of the term is changed to
          54             :     !        a "+".
          55             :     !
          56             :     ! The timestep index (t+1) means that the value of var_zt being used is from
          57             :     ! the next timestep, which is being advanced to in solving the d(var_zt)/dt
          58             :     ! equation.
          59             :     !
          60             :     ! This term is discretized as follows:
          61             :     !
          62             :     ! The values of var_zt are found on the thermodynamic levels, as are the
          63             :     ! values of wm_zt (mean vertical velocity on thermodynamic levels).  The
          64             :     ! variable var_zt is interpolated to the intermediate momentum levels.  The
          65             :     ! derivative of the interpolated values is taken over the central
          66             :     ! thermodynamic level.  The derivative is multiplied by wm_zt at the central
          67             :     ! thermodynamic level to get the desired result.
          68             :     !
          69             :     ! -----var_zt---------------------------------------------- t(k+1)
          70             :     !
          71             :     ! =================var_zt(interp)========================== m(k)
          72             :     !
          73             :     ! -----var_zt---------------------d(var_zt)/dz-----wm_zt--- t(k)
          74             :     !
          75             :     ! =================var_zt(interp)========================== m(k-1)
          76             :     !
          77             :     ! -----var_zt---------------------------------------------- t(k-1)
          78             :     !
          79             :     ! The vertical indices t(k+1), m(k), t(k), m(k-1), and t(k-1) correspond
          80             :     ! with altitudes zt(k+1), zm(k), zt(k), zm(k-1), and zt(k-1), respectively.
          81             :     ! The letter "t" is used for thermodynamic levels and the letter "m" is used
          82             :     ! for momentum levels.
          83             :     !
          84             :     ! invrs_dzt(k) = 1 / ( zm(k) - zm(k-1) )
          85             :     !
          86             :     !
          87             :     ! Special discretization for upper boundary level:
          88             :     !
          89             :     ! Method 1:  Constant derivative method (or "one-sided" method).
          90             :     !
          91             :     ! The values of var_zt are found on the thermodynamic levels, as are the
          92             :     ! values of wm_zt (mean vertical velocity on the thermodynamic levels).  The
          93             :     ! variable var_zt is interpolated to momentum level gr%nz-1, based on
          94             :     ! the values of var_zt at thermodynamic levels gr%nz and gr%nz-1.
          95             :     ! However, the variable var_zt cannot be interpolated to momentum level
          96             :     ! gr%nz.  Rather, a linear extension is used to find the value of var_zt
          97             :     ! at momentum level gr%nz, based on the values of var_zt at thermodynamic
          98             :     ! levels gr%nz and gr%nz-1.  The derivative of the extended and
          99             :     ! interpolated values, d(var_zt)/dz, is taken over the central thermodynamic
         100             :     ! level.  Of course, this derivative will be the same as the derivative of
         101             :     ! var_zt between thermodynamic levels gr%nz and gr%nz-1.  The derivative
         102             :     ! is multiplied by wm_zt at the central thermodynamic level to get the
         103             :     ! desired result.
         104             :     !
         105             :     ! For the following diagram, k = gr%nz, which is the uppermost level of
         106             :     ! the model:
         107             :     !
         108             :     ! =================var_zt(extend)========================== m(k)   Boundary
         109             :     !
         110             :     ! -----var_zt---------------------d(var_zt)/dz-----wm_zt--- t(k)
         111             :     !
         112             :     ! =================var_zt(interp)========================== m(k-1)
         113             :     !
         114             :     ! -----var_zt---------------------------------------------- t(k-1)
         115             :     !
         116             :     !
         117             :     ! Method 2:  Zero derivative method:
         118             :     !            the derivative d(var_zt)/dz over the model top is set to 0.
         119             :     !
         120             :     ! This method corresponds with the "zero-flux" boundary condition option
         121             :     ! for eddy diffusion, where d(var_zt)/dz is set to 0 across the upper
         122             :     ! boundary.
         123             :     !
         124             :     ! In order to discretize the upper boundary condition, consider a new level
         125             :     ! outside the model (thermodynamic level gr%nz+1) just above the upper
         126             :     ! boundary level (thermodynamic level gr%nz).  The value of var_zt at the
         127             :     ! level just outside the model is defined to be the same as the value of
         128             :     ! var_zt at thermodynamic level gr%nz.  Therefore, the value of
         129             :     ! d(var_zt)/dz between the level just outside the model and the uppermost
         130             :     ! thermodynamic level is 0, staying consistent with the zero-flux boundary
         131             :     ! condition option for the eddy diffusion portion of the code.  Therefore,
         132             :     ! the value of var_zt at momentum level gr%nz, which is the upper boundary
         133             :     ! of the model, would be the same as the value of var_zt at the uppermost
         134             :     ! thermodynamic level.
         135             :     !
         136             :     ! The values of var_zt are found on the thermodynamic levels, as are the
         137             :     ! values of wm_zt (mean vertical velocity on the thermodynamic levels).  The
         138             :     ! variable var_zt is interpolated to momentum level gr%nz-1, based on
         139             :     ! the values of var_zt at thermodynamic levels gr%nz and gr%nz-1.  The
         140             :     ! value of var_zt at momentum level gr%nz is set equal to the value of
         141             :     ! var_zt at thermodynamic level gr%nz, as described above.  The derivative
         142             :     ! of the set and interpolated values, d(var_zt)/dz, is taken over the
         143             :     ! central thermodynamic level.  The derivative is multiplied by wm_zt at the
         144             :     ! central thermodynamic level to get the desired result.
         145             :     !
         146             :     ! For the following diagram, k = gr%nz, which is the uppermost level of
         147             :     ! the model:
         148             :     !
         149             :     ! --[var_zt(k+1) = var_zt(k)]----(level outside model)----- t(k+1)
         150             :     !
         151             :     ! ==[var_zt(top) = var_zt(k)]===[d(var_zt)/dz|_(top) = 0]== m(k)   Boundary
         152             :     !
         153             :     ! -----var_zt(k)------------------d(var_zt)/dz-----wm_zt--- t(k)
         154             :     !
         155             :     ! =================var_zt(interp)========================== m(k-1)
         156             :     !
         157             :     ! -----var_zt(k-1)----------------------------------------- t(k-1)
         158             :     !
         159             :     ! where (top) stands for the grid index of momentum level k = gr%nz, which
         160             :     ! is the upper boundary of the model.
         161             :     !
         162             :     ! This method of boundary discretization is also similar to the method
         163             :     ! currently employed at the lower boundary for most thermodynamic-level
         164             :     ! variables.  Since thermodynamic level k = 1 is below the model bottom,
         165             :     ! mean advection is not applied.  Thus, thermodynamic level k = 2 becomes
         166             :     ! the lower boundary level.  Now, the mean advection term at thermodynamic
         167             :     ! level 2 takes into account var_zt from levels 1, 2, and 3.  However, in
         168             :     ! most cases, the value of var_zt(1) is set equal to var_zt(2) after the
         169             :     ! matrix of equations has been solved.  Therefore, the derivative,
         170             :     ! d(var_zt)/dz, over the model bottom (momentum level k = 1) becomes 0.
         171             :     ! Thus, the method of setting d(var_zt)/dz to 0 over the model top keeps
         172             :     ! the way the upper and lower boundaries are handled consistent with each
         173             :     ! other.
         174             :
         175             :     ! References:
         176             :     !   None
         177             :     !-----------------------------------------------------------------------
         178             :
         179             :     use grid_class, only: &
         180             :         grid ! Type
         181             :
         182             :     use constants_clubb, only: &
         183             :         one,  & ! Constant(s)
         184             :         zero
         185             :
         186             :     use clubb_precision, only: &
         187             :         core_rknd ! Variable(s)
         188             :
         189             :     implicit none
         190             :
         191             :     ! Constant parameters
         192             :     integer, parameter :: &
         193             :       kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
         194             :       k_tdiag   = 2, & ! Thermodynamic main diagonal index.
         195             :       km1_tdiag = 3    ! Thermodynamic subdiagonal index.
         196             :
         197             :     integer, parameter :: &
         198             :       t_above = 1, & ! Index for upper thermodynamic level grid weight.
         199             :       t_below = 2    ! Index for lower thermodynamic level grid weight.
         200             :
         201             :     ! -------------------------- Input Variables --------------------------
         202             :     integer, intent(in) :: &
         203             :       nz, &
         204             :       ngrdcol
         205             :
         206             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         207             :       wm_zt,     & ! wm_zt                        [m/s]
         208             :       invrs_dzt, & ! Inverse of grid spacing      [1/m]
         209             :       invrs_dzm    ! Inverse of grid spacing      [1/m]
         210             :
         211             :     real( kind = core_rknd ), dimension(ngrdcol,nz,t_above:t_below), intent(in) :: &
         212             :       weights_zt2zm
         213             :
         214             :     logical, intent(in) :: &
         215             :       l_upwind_xm_ma    ! This flag determines whether we want to use an upwind
         216             :                         ! differencing approximation rather than a centered
         217             :                         ! differencing for turbulent or mean advection terms.
         218             :                         ! It affects rtm, thlm, sclrm, um and vm.
         219             :
         220             :     ! -------------------------- Return Variable --------------------------
         221             :     real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz), intent(out) :: &
         222             :       lhs_ma    ! Mean advection contributions to lhs    [1/s]
         223             :
         224             :     ! -------------------------- Local Variables --------------------------
         225             :
         226             :     integer :: i, k, b    ! Vertical level index
         227             :
         228             :     !-------------------------- Begin Code --------------------------
         229             :
         230             :     !$acc data copyin( wm_zt, invrs_dzt, invrs_dzm, weights_zt2zm ) &
         231             :     !$acc      copyout( lhs_ma )
         232             :
         233             :     ! Set lower boundary array to 0
         234             :     !$acc parallel loop gang vector collapse(2) default(present)
         235   149194656 :     do i = 1, ngrdcol
         236   569973456 :       do b = 1, ndiags3
         237   561038400 :         lhs_ma(b,i,1) = 0.0_core_rknd
         238             :       end do
         239             :     end do
         240             :     !$acc end parallel loop
         241             :
         242             :
         243     8935056 :     if ( .not. l_upwind_xm_ma ) then  ! Use centered differencing
         244             :
         245             :       ! Most of the interior model; normal conditions.
         246             :       !$acc parallel loop gang vector collapse(2) default(present)
         247           0 :       do k = 3, nz-1, 1
         248           0 :         do i = 1, ngrdcol
         249             :
         250             :           ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         251           0 :           lhs_ma(kp1_tdiag,i,k) = + wm_zt(i,k) * invrs_dzt(i,k) * weights_zt2zm(i,k,t_above)
         252             :
         253             :           ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         254             :           lhs_ma(k_tdiag,i,k) = + wm_zt(i,k) * invrs_dzt(i,k) * ( weights_zt2zm(i,k,t_below) &
         255           0 :                               - weights_zt2zm(i,k-1,t_above) )
         256             :
         257             :           ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         258           0 :           lhs_ma(km1_tdiag,i,k) = - wm_zt(i,k) * invrs_dzt(i,k) * weights_zt2zm(i,k-1,t_below)
         259             :
         260             :         end do
         261             :       end do ! k = 3, nz-1, 1
         262             :       !$acc end parallel loop
         263             :
         264             :       ! Upper Boundary
         265             :
         266             :       ! Special discretization for zero derivative method, where the
         267             :       ! derivative d(var_zt)/dz over the model top is set to 0, in order
         268             :       ! to stay consistent with the zero-flux boundary condition option
         269             :       ! in the eddy diffusion code.
         270             :       !$acc parallel loop gang vector default(present)
         271           0 :       do i = 1, ngrdcol
         272             :
         273             :         ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         274           0 :         lhs_ma(kp1_tdiag,i,nz)   = zero
         275             :
         276             :         ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         277             :         lhs_ma(k_tdiag,i,nz) = + wm_zt(i,nz) * invrs_dzt(i,nz) &
         278           0 :                                  * ( one - weights_zt2zm(i,nz-1,t_above) )
         279             :
         280             :         ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         281             :         lhs_ma(km1_tdiag,i,nz) = - wm_zt(i,nz) * invrs_dzt(i,nz) &
         282           0 :                                    * weights_zt2zm(i,nz-1,t_below)
         283             :
         284             :       end do
         285             :       !$acc end parallel loop
         286             :
         287             :       ! Lower Boundary
         288             :
         289             :       ! Special discretization for zero derivative method, where the
         290             :       ! derivative d(var_zt)/dz over the model bottom is set to 0.
         291             :       !$acc parallel loop gang vector default(present)
         292           0 :       do i = 1, ngrdcol
         293             :
         294             :         ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         295           0 :         lhs_ma(kp1_tdiag,i,2) = + wm_zt(i,2) * invrs_dzt(i,2) &
         296           0 :                                   * weights_zt2zm(i,2,t_above)
         297             :
         298             :         ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         299             :         lhs_ma(k_tdiag,i,2) = - wm_zt(i,2) * invrs_dzt(i,2) &
         300           0 :                                 * ( one - weights_zt2zm(i,2,t_below) )
         301             :
         302             :         ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         303           0 :         lhs_ma(km1_tdiag,i,2) = zero
         304             :
         305             :       end do
         306             :       !$acc end parallel loop
         307             :
         308             :     else ! l_upwind_xm_ma == .true.; use "upwind" differencing
         309             :
         310             :       ! Most of the interior model; normal conditions.
         311             :       !$acc parallel loop gang vector collapse(2) default(present)
         312   741609648 :       do k = 3, nz-1, 1
         313 12242896848 :         do i = 1, ngrdcol
         314 12233961792 :           if ( wm_zt(i,k) >= zero ) then  ! Mean wind is in upward direction
         315             :
         316             :              ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         317  5115895680 :              lhs_ma(kp1_tdiag,i,k) = zero
         318             :
         319             :              ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         320  5115895680 :              lhs_ma(k_tdiag,i,k) = + wm_zt(i,k) * invrs_dzm(i,k-1)
         321             :
         322             :              ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         323  5115895680 :              lhs_ma(km1_tdiag,i,k) = - wm_zt(i,k) * invrs_dzm(i,k-1)
         324             :
         325             :           else  ! wm_zt < 0; Mean wind is in downward direction
         326             :
         327             :              ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         328  6385391520 :              lhs_ma(kp1_tdiag,i,k) = + wm_zt(i,k) * invrs_dzm(i,k)
         329             :
         330             :              ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         331  6385391520 :              lhs_ma(k_tdiag,i,k) = - wm_zt(i,k) * invrs_dzm(i,k)
         332             :
         333             :              ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         334  6385391520 :              lhs_ma(km1_tdiag,i,k) = zero
         335             :
         336             :           endif ! wm_zt > 0
         337             :
         338             :         end do
         339             :       end do ! k = 3, nz-1, 1
         340             :       !$acc end parallel loop
         341             :
         342             :       ! Upper Boundary
         343             :       !$acc parallel loop gang vector default(present)
         344   149194656 :       do i = 1, ngrdcol
         345   149194656 :         if ( wm_zt(i,nz) >= zero ) then  ! Mean wind is in upward direction
         346             :
         347             :           ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         348    65298330 :           lhs_ma(kp1_tdiag,i,nz) = zero
         349             :
         350             :           ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         351    65298330 :           lhs_ma(k_tdiag,i,nz) = + wm_zt(i,nz) * invrs_dzm(i,nz-1)
         352             :
         353             :           ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         354    65298330 :           lhs_ma(km1_tdiag,i,nz) = - wm_zt(i,nz) * invrs_dzm(i,nz-1)
         355             :
         356             :         else  ! wm_zt < 0; Mean wind is in downward direction
         357             :
         358             :           ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         359    74961270 :           lhs_ma(kp1_tdiag,i,nz) = zero
         360             :
         361             :           ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         362    74961270 :           lhs_ma(k_tdiag,i,nz) = zero
         363             :
         364             :           ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         365    74961270 :           lhs_ma(km1_tdiag,i,nz) = zero
         366             :
         367             :         end if ! wm_zt > 0
         368             :       end do
         369             :       !$acc end parallel loop
         370             :
         371             :       ! Lower Boundary
         372             :       !$acc parallel loop gang vector default(present)
         373   149194656 :       do i = 1, ngrdcol
         374   149194656 :         if ( wm_zt(i,2) >= zero ) then  ! Mean wind is in upward direction
         375             :
         376             :            ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         377   140259600 :            lhs_ma(kp1_tdiag,i,2) = zero
         378             :
         379             :            ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         380   140259600 :            lhs_ma(k_tdiag,i,2) = zero
         381             :
         382             :            ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         383   140259600 :            lhs_ma(km1_tdiag,i,2) = zero
         384             :
         385             :         else  ! wm_zt < 0; Mean wind is in downward direction
         386             :
         387             :            ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         388           0 :            lhs_ma(kp1_tdiag,i,2) = + wm_zt(i,2) * invrs_dzm(i,2)
         389             :
         390             :            ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         391           0 :            lhs_ma(k_tdiag,i,2) = - wm_zt(i,2) * invrs_dzm(i,2)
         392             :
         393             :            ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         394           0 :            lhs_ma(km1_tdiag,i,2) = zero
         395             :
         396             :         endif ! wm_zt > 0
         397             :       end do
         398             :       !$acc end parallel loop
         399             :
         400             :     endif ! l_upwind_xm_ma
         401             :
         402             :     !$acc end data
         403             :
         404     8935056 :     return
         405             :
         406             :   end subroutine term_ma_zt_lhs
         407             :
         408             :   !=============================================================================
         409    26805168 :   subroutine term_ma_zm_lhs( nz, ngrdcol, wm_zm, &
         410    26805168 :                                   invrs_dzm, weights_zm2zt, &
         411    26805168 :                                   lhs_ma )
         412             :
         413             :     ! Description:
         414             :     ! Mean advection of var_zm:  implicit portion of the code.
         415             :     !
         416             :     ! The variable "var_zm" stands for a variable that is located at momentum
         417             :     ! grid levels.
         418             :     !
         419             :     ! The d(var_zm)/dt equation contains a mean advection term:
         420             :     !
         421             :     ! - w * d(var_zm)/dz.
         422             :     !
         423             :     ! This term is solved for completely implicitly, such that:
         424             :     !
         425             :     ! - w * d( var_zm(t+1) )/dz.
         426             :     !
         427             :     ! Note:  When the term is brought over to the left-hand side, the sign
         428             :     !        is reversed and the leading "-" in front of the term is changed to
         429             :     !        a "+".
         430             :     !
         431             :     ! The timestep index (t+1) means that the value of var_zm being used is from
         432             :     ! the next timestep, which is being advanced to in solving the d(var_zm)/dt
         433             :     ! equation.
         434             :     !
         435             :     ! This term is discretized as follows:
         436             :     !
         437             :     ! The values of var_zm are found on the momentum levels, as are the values
         438             :     ! of wm_zm (mean vertical velocity on momentum levels).  The variable var_zm
         439             :     ! is interpolated to the intermediate thermodynamic levels.  The derivative
         440             :     ! of the interpolated values is taken over the central momentum level.  The
         441             :     ! derivative is multiplied by wm_zm at the central momentum level to get the
         442             :     ! desired result.
         443             :     !
         444             :     ! =====var_zm============================================== m(k+1)
         445             :     !
         446             :     ! -----------------var_zm(interp)-------------------------- t(k+1)
         447             :     !
         448             :     ! =====var_zm=====================d(var_zm)/dz=====wm_zm=== m(k)
         449             :     !
         450             :     ! -----------------var_zm(interp)-------------------------- t(k)
         451             :     !
         452             :     ! =====var_zm============================================== m(k-1)
         453             :     !
         454             :     ! The vertical indices m(k+1), t(k+1), m(k), t(k), and m(k-1) correspond
         455             :     ! with altitudes zm(k+1), zt(k+1), zm(k), zt(k), and zm(k-1), respectively.
         456             :     ! The letter "t" is used for thermodynamic levels and the letter "m" is used
         457             :     ! for momentum levels.
         458             :     !
         459             :     ! invrs_dzm(k) = 1 / ( zt(k+1) - zt(k) )
         460             :
         461             :     ! References:
         462             :     !-----------------------------------------------------------------------
         463             :
         464             :     use grid_class, only: &
         465             :         grid ! Type
         466             :
         467             :     use constants_clubb, only: &
         468             :         zero ! Constant(s)
         469             :
         470             :     use clubb_precision, only: &
         471             :         core_rknd ! Variable(s)
         472             :
         473             :     implicit none
         474             :
         475             :     ! Constant parameters
         476             :     integer, parameter :: &
         477             :       kp1_mdiag = 1, & ! Momentum superdiagonal index.
         478             :       k_mdiag   = 2, & ! Momentum main diagonal index.
         479             :       km1_mdiag = 3    ! Momentum subdiagonal index.
         480             :
         481             :     integer, parameter :: &
         482             :       m_above = 1, & ! Index for upper momentum level grid weight.
         483             :       m_below = 2    ! Index for lower momentum level grid weight.
         484             :
         485             :     ! -------------------------- Input Variables --------------------------
         486             :     integer, intent(in) :: &
         487             :       nz, &
         488             :       ngrdcol
         489             :
         490             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         491             :       wm_zm,     & ! wm_zm                        [m/s]
         492             :       invrs_dzm    ! Inverse of grid spacing      [1/m]
         493             :
         494             :     real( kind = core_rknd ), dimension(ngrdcol,nz,m_above:m_below), intent(in) :: &
         495             :       weights_zm2zt
         496             :
         497             :     ! -------------------------- Return Variable --------------------------
         498             :     real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz), intent(out) :: &
         499             :       lhs_ma    ! Mean advection contributions to lhs  [1/s]
         500             :
         501             :     ! -------------------------- Local Variables
         502             :     integer :: i, k, b    ! Vertical level index
         503             :
         504             :     ! -------------------------- Begin Code --------------------------
         505             :
         506             :     !$acc data copyin( wm_zm, invrs_dzm, weights_zm2zt ) &
         507             :     !$acc      copyout( lhs_ma )
         508             :
         509             :     ! Set lower boundary array to 0
         510             :     !$acc parallel loop gang vector collapse(2) default(present)
         511   447583968 :     do i = 1, ngrdcol
         512  1709920368 :       do b = 1, ndiags3
         513  1683115200 :         lhs_ma(b,i,1) = zero
         514             :       end do
         515             :     end do
         516             :     !$acc end parallel loop
         517             :
         518             :     ! Most of the interior model; normal conditions.
         519             :     !$acc parallel loop gang vector collapse(2) default(present)
         520  2251634112 :     do k = 2, nz-1, 1
         521 37176274512 :       do i = 1, ngrdcol
         522             :
         523             :         ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
         524 34924640400 :         lhs_ma(kp1_mdiag,i,k) = + wm_zm(i,k) * invrs_dzm(i,k) * weights_zm2zt(i,k+1,m_above)
         525             :
         526             :         ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
         527             :         lhs_ma(k_mdiag,i,k) = + wm_zm(i,k) * invrs_dzm(i,k) * ( weights_zm2zt(i,k+1,m_below) &
         528 34924640400 :                                        - weights_zm2zt(i,k,m_above) )
         529             :
         530             :         ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
         531 37149469344 :         lhs_ma(km1_mdiag,i,k) = - wm_zm(i,k) * invrs_dzm(i,k) * weights_zm2zt(i,k,m_below)
         532             :
         533             :       end do
         534             :     end do ! k = 2, nz-1, 1
         535             :     !$acc end parallel loop
         536             :
         537             :     ! Set upper boundary array to 0
         538             :     !$acc parallel loop gang vector collapse(2) default(present)
         539   447583968 :     do i = 1, ngrdcol
         540  1709920368 :       do b = 1, ndiags3
         541  1683115200 :         lhs_ma(b,i,nz) = zero
         542             :       end do
         543             :     end do
         544             :     !$acc end parallel loop
         545             :
         546             :     !$acc end data
         547             :
         548    26805168 :     return
         549             :
         550             :   end subroutine term_ma_zm_lhs
         551             :
         552             :   !=============================================================================
         553             :
         554             : end module mean_adv

Generated by: LCOV version 1.14