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 17:57:11 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  5107411116 :              lhs_ma(kp1_tdiag,i,k) = zero
     318             : 
     319             :              ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
     320  5107411116 :              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  5107411116 :              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  6393876084 :              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  6393876084 :              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  6393876084 :              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    65343546 :           lhs_ma(kp1_tdiag,i,nz) = zero
     349             : 
     350             :           ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
     351    65343546 :           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    65343546 :           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    74916054 :           lhs_ma(kp1_tdiag,i,nz) = zero
     360             : 
     361             :           ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
     362    74916054 :           lhs_ma(k_tdiag,i,nz) = zero
     363             : 
     364             :           ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
     365    74916054 :           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