LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - diffusion.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 74 | 105 | 70.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 3 | 66.7 %

              Line data    Source code
           1             : !-----------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module diffusion
           5             :
           6             :   ! Description:
           7             :   ! Module diffusion computes the eddy diffusion terms for all of the
           8             :   ! time-tendency (prognostic) equations in the CLUBB parameterization.  Most of
           9             :   ! the eddy diffusion terms are solved for completely implicitly, and therefore
          10             :   ! become part of the left-hand side of their respective equations.  However,
          11             :   ! wp2 and wp3 have an option to use a Crank-Nicholson eddy diffusion scheme,
          12             :   ! which has both implicit and explicit components.
          13             :   !
          14             :   ! Function diffusion_zt_lhs handles the eddy diffusion terms for the variables
          15             :   ! located at thermodynamic grid levels.  These variables are:  wp3 and all
          16             :   ! hydrometeor species.  The variables um and vm also use the Crank-Nicholson
          17             :   ! eddy-diffusion scheme for their turbulent advection term.
          18             :   !
          19             :   ! Function diffusion_zm_lhs handles the eddy diffusion terms for the variables
          20             :   ! located at momentum grid levels.  The variables are:  wprtp, wpthlp, wp2,
          21             :   ! rtp2, thlp2, rtpthlp, up2, vp2, wpsclrp, sclrprtp, sclrpthlp, and sclrp2.
          22             :
          23             :   implicit none
          24             :
          25             :   private ! Default Scope
          26             :
          27             :   public :: diffusion_zt_lhs, &
          28             :             diffusion_cloud_frac_zt_lhs, &
          29             :             diffusion_zm_lhs
          30             :
          31             :   contains
          32             :
          33             :   !=============================================================================
          34    17870112 :   subroutine diffusion_zt_lhs( nz, ngrdcol, gr, K_zm, K_zt, nu,  & ! In
          35    17870112 :                                     invrs_rho_ds_zt, rho_ds_zm,       & ! In
          36    17870112 :                                     lhs )                               ! Out
          37             :
          38             :     ! Description:
          39             :     ! Vertical eddy diffusion 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 an eddy diffusion term:
          45             :     !
          46             :     ! + d [ ( K_zm + nu ) * d(var_zt)/dz ] / dz.
          47             :     !
          48             :     ! This term is usually solved for completely implicitly, such that:
          49             :     !
          50             :     ! + d [ ( K_zm + nu ) * d( var_zt(t+1) )/dz ] / dz.
          51             :     !
          52             :     ! However, when a Crank-Nicholson scheme is used, the eddy diffusion term
          53             :     ! has both implicit and explicit components, such that:
          54             :     !
          55             :     ! + (1/2) * d [ ( K_zm + nu ) * d( var_zt(t+1) )/dz ] / dz
          56             :     !    + (1/2) * d [ ( K_zm + nu ) * d( var_zt(t) )/dz ] / dz;
          57             :     !
          58             :     ! for which the implicit component is:
          59             :     !
          60             :     ! + (1/2) * d [ ( K_zm + nu ) * d( var_zt(t+1) )/dz ] / dz.
          61             :     !
          62             :     ! Note:  When the implicit term is brought over to the left-hand side,
          63             :     !        the  sign is reversed and the leading "+" in front of the term
          64             :     !        is changed to a "-".
          65             :     !
          66             :     ! Timestep index (t) stands for the index of the current timestep, while
          67             :     ! timestep index (t+1) stands for the index of the next timestep, which is
          68             :     ! being advanced to in solving the d(var_zt)/dt equation.
          69             :     !
          70             :     ! The implicit portion of this term is discretized as follows:
          71             :     !
          72             :     ! The values of var_zt are found on the thermodynamic levels, while the
          73             :     ! values of K_zm are found on the momentum levels.  The derivatives (d/dz)
          74             :     ! of var_zt are taken over the intermediate momentum levels.  At the
          75             :     ! intermediate momentum levels, d(var_zt)/dz is multiplied by ( K_zm + nu ).
          76             :     ! Then, the derivative of the whole mathematical expression is taken over
          77             :     ! the central thermodynamic level, which yields the desired result.
          78             :     !
          79             :     ! --var_zt------------------------------------------------- t(k+1)
          80             :     !
          81             :     ! ==========d(var_zt)/dz==(K_zm+nu)======================== m(k)
          82             :     !
          83             :     ! --var_zt-------------------d[(K_zm+nu)*d(var_zt)/dz]/dz-- t(k)
          84             :     !
          85             :     ! ==========d(var_zt)/dz==(K_zm+nu)======================== m(k-1)
          86             :     !
          87             :     ! --var_zt------------------------------------------------- t(k-1)
          88             :     !
          89             :     ! The vertical indices t(k+1), m(k), t(k), m(k-1), and t(k-1) correspond
          90             :     ! with altitudes zt(k+1), zm(k), zt(k), zm(k-1), and zt(k-1), respectively.
          91             :     ! The letter "t" is used for thermodynamic levels and the letter "m" is used
          92             :     ! for momentum levels.
          93             :     !
          94             :     ! invrs_dzt(k)   = 1 / ( zm(k) - zm(k-1) )
          95             :     ! invrs_dzm(k)   = 1 / ( zt(k+1) - zt(k) )
          96             :     ! invrs_dzm(k-1) = 1 / ( zt(k) - zt(k-1) )
          97             :     !
          98             :     ! Note:  This function only computes the general implicit form:
          99             :     !        + d [ ( K_zm + nu ) * d( var_zt(t+1) )/dz ] / dz.
         100             :     !        For a Crank-Nicholson scheme, the left-hand side result of this
         101             :     !        function will have to be multiplied by (1/2).  For a
         102             :     !        Crank-Nicholson scheme, the right-hand side (explicit) component
         103             :     !        needs to be computed by multiplying the left-hand side results by
         104             :     !        (1/2), reversing the sign on each left-hand side element, and then
         105             :     !        multiplying each element by the appropriate var_zt(t) value from
         106             :     !        the appropriate vertical level.
         107             :     !
         108             :     !
         109             :     ! Boundary Conditions:
         110             :     !
         111             :     ! 1) Zero-flux boundary conditions.
         112             :     !    This function is set up to use zero-flux boundary conditions at both
         113             :     !    the lower boundary level and the upper boundary level.  The flux, F,
         114             :     !    is the amount of var_zt flowing normal through the boundary per unit
         115             :     !    time per unit surface area.  The derivative of the flux effects the
         116             :     !    time-tendency of var_zt, such that:
         117             :     !
         118             :     !    d(var_zt)/dt = -dF/dz.
         119             :     !
         120             :     !    For the 2nd-order eddy-diffusion term, +d[(K_zm+nu)*d(var_zt)/dz]/dz,
         121             :     !    the flux is:
         122             :     !
         123             :     !    F = -(K_zm+nu)*d(var_zt)/dz.
         124             :     !
         125             :     !    In order to have zero-flux boundary conditions, the derivative of
         126             :     !    var_zt, d(var_zt)/dz, needs to equal 0 at both the lower boundary and
         127             :     !    the upper boundary.
         128             :     !
         129             :     !    In order to discretize the lower boundary condition, consider a new
         130             :     !    level outside the model (thermodynamic level 0) just below the lower
         131             :     !    boundary level (thermodynamic level 1).  The value of var_zt at the
         132             :     !    level just outside the model is defined to be the same as the value of
         133             :     !    var_zt at the lower boundary level.  Therefore, the value of
         134             :     !    d(var_zt)/dz between the level just outside the model and the lower
         135             :     !    boundary level is 0, satisfying the zero-flux boundary condition.  The
         136             :     !    other value for d(var_zt)/dz (between thermodynamic level 2 and
         137             :     !    thermodynamic level 1) is taken over the intermediate momentum level
         138             :     !    (momentum level 1), where it is multiplied by the factor
         139             :     !    ( K_zm(1) + nu ).  Then, the derivative of the whole expression is
         140             :     !    taken over the central thermodynamic level.
         141             :     !
         142             :     !    -var_zt(2)-------------------------------------------- t(2)
         143             :     !
         144             :     !    ==========d(var_zt)/dz==(K_zm(1)+nu)================== m(1)
         145             :     !
         146             :     !    -var_zt(1)---------------d[(K_zm+nu)*d(var_zt)/dz]/dz- t(1) Boundary
         147             :     !
         148             :     !              [d(var_zt)/dz = 0]
         149             :     !
         150             :     !    -[var_zt(0) = var_zt(1)]-----(level outside model)---- t(0)
         151             :     !
         152             :     !    The result is dependent only on values of K_zm found at momentum
         153             :     !    level 1 and values of var_zt found at thermodynamic levels 1 and 2.
         154             :     !    Thus, it only affects 2 diagonals on the left-hand side matrix.
         155             :     !
         156             :     !    The same method can be used to discretize the upper boundary by
         157             :     !    considering a new level outside the model just above the upper boundary
         158             :     !    level.
         159             :     !
         160             :     ! 2) Fixed-point boundary conditions.
         161             :     !    Many equations in the model use fixed-point boundary conditions rather
         162             :     !    than zero-flux boundary conditions.  This means that the value of
         163             :     !    var_zt stays the same over the course of the timestep at the lower
         164             :     !    boundary, as well as at the upper boundary.
         165             :     !
         166             :     !    In order to discretize the boundary conditions for equations requiring
         167             :     !    fixed-point boundary conditions, either:
         168             :     !    a) in the parent subroutine or function (that calls this function),
         169             :     !       loop over all vertical levels from the second-lowest to the
         170             :     !       second-highest, ignoring the boundary levels.  Then set the values
         171             :     !       at the boundary levels in the parent subroutine; or
         172             :     !    b) in the parent subroutine or function, loop over all vertical levels
         173             :     !       and then overwrite the results at the boundary levels.
         174             :     !
         175             :     !    Either way, at the boundary levels, an array with a value of 1 at the
         176             :     !    main diagonal on the left-hand side and with values of 0 at all other
         177             :     !    diagonals on the left-hand side will preserve the right-hand side value
         178             :     !    at that level, thus satisfying the fixed-point boundary conditions.
         179             :     !
         180             :     !
         181             :     ! Conservation Properties:
         182             :     !
         183             :     ! When zero-flux boundary conditions are used, this technique of
         184             :     ! discretizing the eddy diffusion term leads to conservative differencing.
         185             :     ! When conservative differencing is in place, the column totals for each
         186             :     ! column in the left-hand side matrix (for the eddy diffusion term) should
         187             :     ! be equal to 0.  This ensures that the total amount of the quantity var_zt
         188             :     ! over the entire vertical domain is being conserved, meaning that nothing
         189             :     ! is lost due to diffusional effects.
         190             :     !
         191             :     ! To see that this conservation law is satisfied, compute the eddy diffusion
         192             :     ! of var_zt and integrate vertically.  In discretized matrix notation (where
         193             :     ! "i" stands for the matrix column and "j" stands for the matrix row):
         194             :     !
         195             :     !  0 = Sum_j Sum_i ( 1/invrs_dzt )_i
         196             :     !                     ( invrs_dzt * ((K_zm+nu)*invrs_dzm) )_ij (var_zt)_j.
         197             :     !
         198             :     ! The left-hand side matrix, ( invrs_dzt * ((K_zm+nu)*invrs_dzm) )_ij, is
         199             :     ! partially written below.  The sum over i in the above equation removes
         200             :     ! invrs_dzt everywhere from the matrix below.  The sum over j leaves the
         201             :     ! column totals that are desired.
         202             :     !
         203             :     ! Left-hand side matrix contributions from eddy diffusion term; first four
         204             :     ! vertical levels:
         205             :     !
         206             :     !     -------------------------------------------------------------------------->
         207             :     !k=1 | +invrs_dzt(k)          -invrs_dzt(k)                          0
         208             :     !    |   *(K_zm(k)+nu)          *(K_zm(k)+nu)
         209             :     !    |     *invrs_dzm(k)          *invrs_dzm(k)
         210             :     !    |
         211             :     !k=2 | -invrs_dzt(k)          +invrs_dzt(k)              -invrs_dzt(k)
         212             :     !    |   *(K_zm(k-1)+nu)        *[ (K_zm(k)+nu)            *(K_zm(k)+nu)
         213             :     !    |     *invrs_dzm(k-1)          *invrs_dzm(k)            *invrs_dzm(k)
         214             :     !    |                            +(K_zm(k-1)+nu)
         215             :     !    |                              *invrs_dzm(k-1) ]
         216             :     !    |
         217             :     !k=3 |         0              -invrs_dzt(k)              +invrs_dzt(k)
         218             :     !    |                          *(K_zm(k-1)+nu)            *[ (K_zm(k)+nu)
         219             :     !    |                            *invrs_dzm(k-1)              *invrs_dzm(k)
         220             :     !    |                                                       +(K_zm(k-1)+nu)
         221             :     !    |                                                         *invrs_dzm(k-1) ]
         222             :     !    |
         223             :     !k=4 |         0                        0                -invrs_dzt(k)
         224             :     !    |                                                     *(K_zm(k-1)+nu)
         225             :     !    |                                                       *invrs_dzm(k-1)
         226             :     !   \ /
         227             :     !
         228             :     ! Note:  The superdiagonal term from level 3 and both the main diagonal and
         229             :     !        superdiagonal terms from level 4 are not shown on this diagram.
         230             :     !
         231             :     ! Note:  The matrix shown is a tridiagonal matrix.  For a band diagonal
         232             :     !        matrix (with 5 diagonals), there would be an extra row between each
         233             :     !        of the rows shown and an extra column between each of the columns
         234             :     !        shown.  However, for the purposes of the var_zt eddy diffusion
         235             :     !        term, those extra row and column values are all 0, and the
         236             :     !        conservation properties of the matrix aren't effected.
         237             :     !
         238             :     ! If fixed-point boundary conditions are used, the matrix entries at
         239             :     ! level 1 (k=1) read:  1   0   0; which means that conservative differencing
         240             :     ! is not in play.  The total amount of var_zt over the entire vertical
         241             :     ! domain is not being conserved, as amounts of var_zt may be fluxed out
         242             :     ! through the upper boundary or lower boundary through the effects of
         243             :     ! diffusion.
         244             :     !
         245             :     ! Brian Griffin.  April 26, 2008.
         246             :
         247             :     ! References:
         248             :     ! None
         249             :     !-----------------------------------------------------------------------
         250             :
         251             :     use grid_class, only: &
         252             :         grid, &  ! Type
         253             :         ddzm     ! Procedure
         254             :
         255             :     use constants_clubb, only: &
         256             :         zero    ! Constant(s)
         257             :
         258             :     use clubb_precision, only: &
         259             :         core_rknd ! Variable(s)
         260             :
         261             :     use model_flags, only: &
         262             :         l_upwind_Kh_dp_term
         263             :
         264             :     implicit none
         265             :
         266             :     ! Constant parameters
         267             :     integer, parameter :: &
         268             :       kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
         269             :       k_tdiag   = 2, & ! Thermodynamic main diagonal index.
         270             :       km1_tdiag = 3    ! Thermodynamic subdiagonal index.
         271             :
         272             :     !------------------------ Input Variables ------------------------
         273             :     integer, intent(in) :: &
         274             :       nz, &
         275             :       ngrdcol
         276             :
         277             :     type (grid), target, intent(in) :: gr
         278             :
         279             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         280             :       K_zm,            & ! Coef. of eddy diffusivity at momentum levels   [m^2/s]
         281             :       K_zt,            & ! Coef. of eddy diffusivity at thermo. levels    [m^2/s]
         282             :       rho_ds_zm,       & ! Dry statis density on momentum levels          [kg]
         283             :       invrs_rho_ds_zt    ! Inverse dry statis density on thermo. levels   [1/kg]
         284             :
         285             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
         286             :       nu                ! Background constant coef. of eddy diffusivity  [m^2/s]
         287             :
         288             :     !------------------------ Return Variable ------------------------
         289             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(out) :: &
         290             :       lhs     ! LHS coefficient of diffusion term    [1/s]
         291             :
         292             :     !------------------------ Local Variable ------------------------
         293             :     integer :: i,k    ! Vertical level index
         294             :
         295             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz) :: &
         296    17870112 :       lhs_upwind    ! LHS coefficient diffusion term due to upwind  [1/s]
         297             :
         298             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         299    17870112 :       drhoKdz_zt, &
         300    35740224 :       K_zm_nu, &
         301    17870112 :       rho_K_zm_nu, &
         302    17870112 :       ddzm_rho_K_zm_nu
         303             :
         304             :     !------------------------ Begin code ------------------------
         305             :
         306             :     !$acc data copyin( gr, gr%invrs_dzm, gr%invrs_dzt, &
         307             :     !$acc              K_zm, K_zt, rho_ds_zm, invrs_rho_ds_zt, nu ) &
         308             :     !$acc     copyout( lhs ) &
         309             :     !$acc      create( lhs_upwind, drhoKdz_zt, K_zm_nu, rho_K_zm_nu, ddzm_rho_K_zm_nu )
         310             :
         311             :     !$acc parallel loop gang vector collapse(2) default(present)
         312  1536829632 :     do k = 1, nz
         313 25380961632 :       do i = 1, ngrdcol
         314 25363091520 :         K_zm_nu(i,k) = K_zm(i,k) + nu(i)
         315             :       end do
         316             :     end do
         317             :     !$acc end parallel loop
         318             :
         319             :     if ( l_upwind_Kh_dp_term ) then
         320             :
         321             :       !$acc parallel loop gang vector collapse(2) default(present)
         322             :       do k = 1, nz
         323             :         do i = 1, ngrdcol
         324             :           rho_K_zm_nu(i,k) = rho_ds_zm(i,k) * K_zm_nu(i,k)
         325             :         end do
         326             :       end do
         327             :       !$acc end parallel loop
         328             :
         329             :       ! calculate the dKh_zt/dz
         330             :       ddzm_rho_K_zm_nu = ddzm( nz, ngrdcol, gr, rho_K_zm_nu )
         331             :
         332             :       !$acc parallel loop gang vector collapse(2) default(present)
         333             :       do k = 1, nz
         334             :         do i = 1, ngrdcol
         335             :           drhoKdz_zt(i,k) = - invrs_rho_ds_zt(i,k) * ddzm_rho_K_zm_nu(i,k)
         336             :         end do
         337             :       end do
         338             :       !$acc end parallel loop
         339             :
         340             :
         341             :       ! extra terms with upwind scheme
         342             :       ! k = 1 (bottom level); lower boundary level
         343             :       !$acc parallel loop gang vector default(present)
         344             :       do i = 1, ngrdcol
         345             :         lhs_upwind(kp1_tdiag,i,1) = zero
         346             :         lhs_upwind(k_tdiag,i,1)   = zero
         347             :         lhs_upwind(km1_tdiag,i,1) = zero
         348             :       end do
         349             :       !$acc end parallel loop
         350             :
         351             :       ! extra terms with upwind scheme
         352             :       ! k = 2 (bottom level); lower boundary level
         353             :       !$acc parallel loop gang vector default(present)
         354             :       do i = 1, ngrdcol
         355             :         lhs_upwind(kp1_tdiag,i,2) = + min( drhoKdz_zt(i,2) , zero ) * gr%invrs_dzm(i,2)
         356             :         lhs_upwind(k_tdiag,i,2)   = - min( drhoKdz_zt(i,2) , zero ) * gr%invrs_dzm(i,2)
         357             :         lhs_upwind(km1_tdiag,i,2) = zero
         358             :       end do
         359             :       !$acc end parallel loop
         360             :
         361             :       ! Most of the interior model; normal conditions.
         362             :         !$acc parallel loop gang vector collapse(2) default(present)
         363             :       do k = 3, nz-1, 1
         364             :         do i = 1, ngrdcol
         365             :
         366             :           ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         367             :           lhs_upwind(kp1_tdiag,i,k) = + min( drhoKdz_zt(i,k) , zero ) * gr%invrs_dzm(i,k)
         368             :
         369             :           ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         370             :           lhs_upwind(k_tdiag,i,k)   = - min( drhoKdz_zt(i,k) , zero ) * gr%invrs_dzm(i,k) &
         371             :                                       + max( drhoKdz_zt(i,k) , zero ) * gr%invrs_dzm(i,k-1)
         372             :
         373             :           ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         374             :           lhs_upwind(km1_tdiag,i,k) = - max( drhoKdz_zt(i,k) , zero ) * gr%invrs_dzm(i,k-1)
         375             :
         376             :         end do
         377             :       end do
         378             :       !$acc end parallel loop
         379             :
         380             :       ! k = nz (top level); upper boundary level.
         381             :       ! Only relevant if zero-flux boundary conditions are used.
         382             :       !$acc parallel loop gang vector default(present)
         383             :       do i = 1, ngrdcol
         384             :         lhs_upwind(kp1_tdiag,i,nz) =  zero
         385             :         lhs_upwind(k_tdiag,i,nz)   = + max( drhoKdz_zt(i,nz) , zero ) * gr%invrs_dzm(i,nz-1)
         386             :         lhs_upwind(km1_tdiag,i,nz) = - max( drhoKdz_zt(i,nz) , zero ) * gr%invrs_dzm(i,nz-1)
         387             :       end do
         388             :       !$acc end parallel loop
         389             :
         390             :     end if
         391             :
         392             :     !$acc parallel loop gang vector default(present)
         393   298389312 :     do i = 1, ngrdcol
         394   280519200 :       lhs(kp1_tdiag,i,1) = zero
         395   280519200 :       lhs(k_tdiag,i,1)   = zero
         396   298389312 :       lhs(km1_tdiag,i,1) = zero
         397             :     end do
         398             :     !$acc end parallel loop
         399             :
         400             :     ! k = 2 (bottom level); lower boundary level.
         401             :     ! Only relevant if zero-flux boundary conditions are used.
         402             :     ! These k=2 lines currently do not have any effect on model results.
         403             :     ! These k=2 level of this "lhs" array is not fed into
         404             :     ! the final LHS matrix that will be used to solve for the next timestep.
         405             :
         406             :     if ( .not. l_upwind_Kh_dp_term ) then
         407             :
         408             :       !$acc parallel loop gang vector default(present)
         409   298389312 :       do i = 1, ngrdcol
         410             :         ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         411   561038400 :         lhs(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
         412   561038400 :                                * ( K_zm(i,2) + nu(i) ) * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
         413             :
         414             :         ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         415           0 :         lhs(k_tdiag,i,2)   = + gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
         416   280519200 :                                * ( K_zm(i,2) + nu(i) ) * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
         417             :
         418             :         ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         419   298389312 :         lhs(km1_tdiag,i,2) = zero
         420             :       end do
         421             :       !$acc end parallel loop
         422             :
         423             :     else
         424             :
         425             :       !$acc parallel loop gang vector default(present)
         426             :       do i = 1, ngrdcol
         427             :         ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         428             :         lhs(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) * ( K_zt(i,2) + nu(i) ) * gr%invrs_dzm(i,2) &
         429             :                              + lhs_upwind(kp1_tdiag,i,2)
         430             :
         431             :         ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         432             :         lhs(k_tdiag,i,2)   = + gr%invrs_dzt(i,2) * ( K_zt(i,2) + nu(i) ) * gr%invrs_dzm(i,2) &
         433             :                              + lhs_upwind(k_tdiag,i,1)
         434             :
         435             :         ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         436             :         lhs(km1_tdiag,i,2) = zero + lhs_upwind(km1_tdiag,i,2)
         437             :       end do
         438             :       !$acc end parallel loop
         439             :
         440             :     end if
         441             :
         442             :     if ( .not. l_upwind_Kh_dp_term ) then
         443             :
         444             :       ! Most of the interior model; normal conditions.
         445             :       !$acc parallel loop gang vector collapse(2) default(present)
         446  1483219296 :       do k = 3, nz-1, 1
         447 24485793696 :         do i = 1, ngrdcol
         448             :
         449             :           ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         450 46005148800 :           lhs(kp1_tdiag,i,k) &
         451           0 :           = - gr%invrs_dzt(i,k) * invrs_rho_ds_zt(i,k) &
         452 69007723200 :                            * K_zm_nu(i,k) * rho_ds_zm(i,k) * gr%invrs_dzm(i,k)
         453             :
         454             :           ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         455             :           lhs(k_tdiag,i,k) &
         456           0 :           = + gr%invrs_dzt(i,k) * invrs_rho_ds_zt(i,k) &
         457           0 :                            * ( K_zm_nu(i,k) * rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
         458 23002574400 :                                + K_zm_nu(i,k-1) * rho_ds_zm(i,k-1) * gr%invrs_dzm(i,k-1) )
         459             :
         460             :           ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         461             :           lhs(km1_tdiag,i,k) &
         462           0 :           = - gr%invrs_dzt(i,k) * invrs_rho_ds_zt(i,k) &
         463 24467923584 :                            * K_zm_nu(i,k-1) * rho_ds_zm(i,k-1) * gr%invrs_dzm(i,k-1)
         464             :         end do
         465             :       end do ! k = 2, nz-1, 1
         466             :       !$acc end parallel loop
         467             :
         468             :     else
         469             :
         470             :       !$acc parallel loop gang vector collapse(2) default(present)
         471             :       do k = 2, nz-1, 1
         472             :         do i = 1, ngrdcol
         473             :           ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         474             :           lhs(kp1_tdiag,i,k) &
         475             :           = - gr%invrs_dzt(i,k) * ( K_zt(i,k) + nu(i) ) * gr%invrs_dzm(i,k) &
         476             :             + lhs_upwind(kp1_tdiag,i,k)
         477             :
         478             :           ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         479             :           lhs(k_tdiag,i,k) &
         480             :           = + gr%invrs_dzt(i,k) * ( K_zt(i,k) + nu(i) ) * ( gr%invrs_dzm(i,k) + gr%invrs_dzm(i,k-1) ) &
         481             :             + lhs_upwind(k_tdiag,i,k)
         482             :
         483             :           ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         484             :           lhs(km1_tdiag,i,k) &
         485             :           = - gr%invrs_dzt(i,k) * ( K_zt(i,k) + nu(i) ) * gr%invrs_dzm(i,k-1) &
         486             :             + lhs_upwind(km1_tdiag,i,k)
         487             :         end do
         488             :       end do ! k = 2, nz-1, 1
         489             :       !$acc end parallel loop
         490             :
         491             :     end if
         492             :
         493             :     ! k = nz (top level); upper boundary level.
         494             :     ! Only relevant if zero-flux boundary conditions are used.
         495             :
         496             :     if ( .not. l_upwind_Kh_dp_term ) then
         497             :
         498             :       !$acc parallel loop gang vector default(present)
         499   298389312 :       do i = 1, ngrdcol
         500             :         ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         501   280519200 :         lhs(kp1_tdiag,i,nz) = zero
         502             :
         503             :         ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         504             :         lhs(k_tdiag,i,nz) &
         505           0 :         = + gr%invrs_dzt(i,nz) * invrs_rho_ds_zt(i,nz) &
         506   280519200 :             * ( K_zm(i,nz-1) + nu(i) ) * rho_ds_zm(i,nz-1) * gr%invrs_dzm(i,nz-1)
         507             :
         508             :         ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         509             :         lhs(km1_tdiag,i,nz) &
         510           0 :         = - gr%invrs_dzt(i,nz) * invrs_rho_ds_zt(i,nz) &
         511   298389312 :             * ( K_zm(i,nz-1) + nu(i) ) * rho_ds_zm(i,nz-1) * gr%invrs_dzm(i,nz-1)
         512             :       end do
         513             :       !$acc end parallel loop
         514             :
         515             :     else
         516             :
         517             :       !$acc parallel loop gang vector default(present)
         518             :       do i = 1, ngrdcol
         519             :         ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         520             :         lhs(kp1_tdiag,i,nz) = zero + lhs_upwind(kp1_tdiag,i,nz)
         521             :
         522             :         ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         523             :         lhs(k_tdiag,i,nz) &
         524             :         = + gr%invrs_dzt(i,nz) * ( K_zt(i,nz) + nu(i) ) * gr%invrs_dzm(i,nz-1) &
         525             :           + lhs_upwind(k_tdiag,i,nz)
         526             :
         527             :         ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         528             :         lhs(km1_tdiag,i,nz) &
         529             :         = - gr%invrs_dzt(i,nz) * ( K_zt(i,nz) + nu(i) ) * gr%invrs_dzm(i,nz-1) &
         530             :           + lhs_upwind(km1_tdiag,i,nz)
         531             :       end do
         532             :       !$acc end parallel loop
         533             :
         534             :     end if
         535             :
         536             :     !$acc end data
         537             :
         538    17870112 :     return
         539             :
         540             :   end subroutine diffusion_zt_lhs
         541             :
         542             :   !=============================================================================
         543           0 :   function diffusion_cloud_frac_zt_lhs &
         544             :                 ( gr, K_zm, K_zmm1, cloud_frac_zt, cloud_frac_ztm1, &
         545             :                   cloud_frac_ztp1, cloud_frac_zm, &
         546             :                   cloud_frac_zmm1, &
         547             :                   nu, invrs_dzmm1, invrs_dzm, invrs_dzt, level ) &
         548           0 :   result( lhs )
         549             :
         550             :   ! Description:
         551             :   !   This function adds a weight of cloud fraction to the existing diffusion
         552             :   !   function for number concentration variables (e.g. cloud droplet number
         553             :   !   concentration).  This code should be considered experimental and may
         554             :   !   contain bugs.
         555             :   ! References:
         556             :   !   This algorithm uses equations derived from Guo, et al. 2009.
         557             :   !-----------------------------------------------------------------------------
         558             :
         559             :      use grid_class, only: &
         560             :         grid ! Type
         561             :
         562             :      use clubb_precision, only: &
         563             :          core_rknd ! Variable(s)
         564             :
         565             :       implicit none
         566             :
         567             :     type (grid), target, intent(in) :: gr
         568             :
         569             :     ! External
         570             :       intrinsic :: min
         571             :
         572             :     ! Constant parameters
         573             :       real( kind = core_rknd ), parameter :: &
         574             :       cf_ratio = 10._core_rknd ! Maximum cloud-fraction coefficient applied to Kh_zm
         575             :
         576             :     integer, parameter :: &
         577             :       kp1_tdiag = 1,    & ! Thermodynamic superdiagonal index.
         578             :       k_tdiag   = 2,    & ! Thermodynamic main diagonal index.
         579             :       km1_tdiag = 3       ! Thermodynamic subdiagonal index.
         580             :
         581             :     ! Input Variables
         582             :     real( kind = core_rknd ), intent(in) ::  &
         583             :       K_zm,            & ! Coef. of eddy diffusivity at mom. level (k)   [m^2/s]
         584             :       K_zmm1,          & ! Coef. of eddy diffusivity at mom. level (k-1) [m^2/s]
         585             :       cloud_frac_zt,   & ! Cloud fraction at the thermo. level (k)       [-]
         586             :       cloud_frac_ztm1, & ! Cloud fraction at the thermo. level (k-1)     [-]
         587             :       cloud_frac_ztp1, & ! Cloud fraction at the thermo. level (k+1)     [-]
         588             :       cloud_frac_zm,   & ! Cloud fraction at the momentum level (k)      [-]
         589             :       cloud_frac_zmm1, & ! Cloud fraction at the momentum level (k-1)    [-]
         590             :       invrs_dzt,       & ! Inverse of grid spacing over thermo. lev. (k) [1/m]
         591             :       invrs_dzm,       & ! Inverse of grid spacing over mom. level (k)   [1/m]
         592             :       invrs_dzmm1        ! Inverse of grid spacing over mom. level (k-1) [1/m]
         593             :
         594             :     real( kind = core_rknd ), intent(in) :: &
         595             :       nu                 ! Background constant coef. of eddy diffusivity [m^2/s]
         596             :
         597             :     integer, intent(in) ::  &
         598             :       level     ! Thermodynamic level where calculation occurs.           [-]
         599             :
         600             :     ! Return Variable
         601             :     real( kind = core_rknd ), dimension(3) :: lhs
         602             :
         603             :     ! ---- Begin Code ----
         604             :
         605           0 :     if ( level == 1 ) then
         606             :
         607             :       ! k = 1 (bottom level); lower boundary level.
         608             :       ! Only relevant if zero-flux boundary conditions are used.
         609             :
         610             :       ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         611             : !     lhs(kp1_tdiag) = - invrs_dzt &
         612             : !                        * (K_zm+nu) &
         613             : !                           * ( cloud_frac_zm / cloud_frac_ztp1 ) * invrs_dzm
         614             :       lhs(kp1_tdiag) = - invrs_dzt &
         615             :                          * (K_zm &
         616             :                             * min( cloud_frac_zm / cloud_frac_ztp1, cf_ratio ) &
         617           0 :                             + nu) * invrs_dzm
         618             :
         619             :       ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         620             : !     lhs(k_tdiag)   = + invrs_dzt &
         621             : !                        * (K_zm+nu) &
         622             : !                           * ( cloud_frac_zm / cloud_frac_ztp1 ) * invrs_dzm
         623             :       lhs(k_tdiag)   = + invrs_dzt &
         624             :                          * (K_zm &
         625             :                             * min( cloud_frac_zm / cloud_frac_ztp1, cf_ratio ) &
         626           0 :                             + nu) * invrs_dzm
         627             :
         628             :       ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         629           0 :       lhs(km1_tdiag) = 0.0_core_rknd
         630             :
         631             :
         632           0 :     else if ( level > 1 .and. level < gr%nz ) then
         633             :
         634             :       ! Most of the interior model; normal conditions.
         635             :
         636             :       ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         637             : !     lhs(kp1_tdiag) = - invrs_dzt &
         638             : !                        * (K_zm+nu) &
         639             : !                           * ( cloud_frac_zm / cloud_frac_ztp1 ) * invrs_dzm
         640             : !     lhs(kp1_tdiag) = - invrs_dzt &
         641             : !                        * (K_zm &
         642             : !                           * ( cloud_frac_zm / cloud_frac_ztp1 ) &
         643             : !                           + nu ) * invrs_dzm
         644             :       lhs(kp1_tdiag) = - invrs_dzt &
         645             :                          * (K_zm &
         646             :                             * min( cloud_frac_zm / cloud_frac_ztp1, cf_ratio ) &
         647           0 :                             + nu ) * invrs_dzm
         648             :
         649             :       ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         650             : !     lhs(k_tdiag)   = + invrs_dzt &
         651             : !                        * (  ((K_zm+nu)*cloud_frac_zm)*invrs_dzm &
         652             : !                           + ((K_zmm1+nu)*cloud_frac_zmm1)*invrs_dzmm1 ) &
         653             : !                        / cloud_frac_zt
         654             : !     lhs(k_tdiag)   = + invrs_dzt &
         655             : !                        * ( nu*(invrs_dzm+invrs_dzmm1) + &
         656             : !                               ( ((K_zm*cloud_frac_zm)*invrs_dzm +
         657             : !                                  (K_zmm1*cloud_frac_zmm1)*invrs_dzmm1)&
         658             : !                                / cloud_frac_zt &
         659             : !                               ) &
         660             : !                          )
         661             :       lhs(k_tdiag)   = + invrs_dzt &
         662             :                           * ( nu*(invrs_dzm+invrs_dzmm1) + &
         663             :                               (   K_zm*invrs_dzm* &
         664             :                                      min( cloud_frac_zm / cloud_frac_zt, &
         665             :                                           cf_ratio ) &
         666             :                                 + K_zmm1*invrs_dzmm1* &
         667             :                                      min( cloud_frac_zmm1 / cloud_frac_zt, &
         668             :                                           cf_ratio ) &
         669             :                               ) &
         670           0 :                             )
         671             :
         672             :       ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         673             : !     lhs(km1_tdiag) = - invrs_dzt * (K_zmm1+nu) * &
         674             : !                        ( cloud_frac_zmm1 / cloud_frac_ztm1 ) * invrs_dzmm1
         675             :       lhs(km1_tdiag) = - invrs_dzt &
         676             :                           * (K_zmm1 &
         677             :                              * min( cloud_frac_zmm1 / cloud_frac_ztm1, &
         678             :                                     cf_ratio ) &
         679           0 :                              + nu ) * invrs_dzmm1
         680             :
         681           0 :     else if ( level == gr%nz ) then
         682             :
         683             :       ! k = gr%nz (top level); upper boundary level.
         684             :       ! Only relevant if zero-flux boundary conditions are used.
         685             :
         686             :       ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
         687           0 :       lhs(kp1_tdiag) = 0.0_core_rknd
         688             :
         689             :       ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
         690             : !     lhs(k_tdiag)   = + invrs_dzt &
         691             : !                         *(K_zmm1+nu) &
         692             : !                           *( cloud_frac_zmm1 / cloud_frac_ztm1 ) * invrs_dzmm1
         693             :       lhs(k_tdiag)   = + invrs_dzt &
         694             :                           * (K_zmm1 &
         695             :                              * min( cloud_frac_zmm1 / cloud_frac_ztm1, &
         696             :                                     cf_ratio ) &
         697           0 :                              + nu) * invrs_dzmm1
         698             :
         699             :       ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
         700             : !     lhs(km1_tdiag) = - invrs_dzt * (K_zmm1+nu) * &
         701             : !                        ( cloud_frac_zmm1 / cloud_frac_ztm1 ) * invrs_dzmm1
         702             :       lhs(km1_tdiag) = - invrs_dzt &
         703             :                           * (K_zmm1 &
         704             :                              * min( cloud_frac_zmm1 / cloud_frac_ztm1, &
         705             :                                     cf_ratio ) &
         706           0 :                              + nu) * invrs_dzmm1
         707             :
         708             :     end if
         709             :
         710           0 :     return
         711           0 :   end function diffusion_cloud_frac_zt_lhs
         712             :
         713             :   !=============================================================================
         714    35740224 :   subroutine diffusion_zm_lhs( nz, ngrdcol, gr, K_zt, K_zm, nu,   & ! In
         715    35740224 :                                     invrs_rho_ds_zm, rho_ds_zt, & ! In
         716    35740224 :                                     lhs )                   ! Out
         717             :
         718             :     ! Description:
         719             :     ! Vertical eddy diffusion of var_zm:  implicit portion of the code.
         720             :     !
         721             :     ! The variable "var_zm" stands for a variable that is located at momentum
         722             :     ! grid levels.
         723             :     !
         724             :     ! The d(var_zm)/dt equation contains an eddy diffusion term:
         725             :     !
         726             :     ! + d [ ( K_zt + nu ) * d(var_zm)/dz ] / dz.
         727             :     !
         728             :     ! This term is usually solved for completely implicitly, such that:
         729             :     !
         730             :     ! + d [ ( K_zt + nu ) * d( var_zm(t+1) )/dz ] / dz.
         731             :     !
         732             :     ! However, when a Crank-Nicholson scheme is used, the eddy diffusion term
         733             :     ! has both implicit and explicit components, such that:
         734             :     !
         735             :     ! + (1/2) * d [ ( K_zt + nu ) * d( var_zm(t+1) )/dz ] / dz
         736             :     !    + (1/2) * d [ ( K_zt + nu ) * d( var_zm(t) )/dz ] / dz;
         737             :     !
         738             :     ! for which the implicit component is:
         739             :     !
         740             :     ! + (1/2) * d [ ( K_zt + nu ) * d( var_zm(t+1) )/dz ] / dz.
         741             :     !
         742             :     ! Note:  When the implicit term is brought over to the left-hand side,
         743             :     !        the sign is reversed and the leading "+" in front of the term
         744             :     !        is changed to a "-".
         745             :     !
         746             :     ! Timestep index (t) stands for the index of the current timestep, while
         747             :     ! timestep index (t+1) stands for the index of the next timestep, which is
         748             :     ! being advanced to in solving the d(var_zm)/dt equation.
         749             :     !
         750             :     ! The implicit portion of this term is discretized as follows:
         751             :     !
         752             :     ! The values of var_zm are found on the momentum levels, while the values of
         753             :     ! K_zt are found on the thermodynamic levels.  The derivatives (d/dz) of
         754             :     ! var_zm are taken over the intermediate thermodynamic levels.  At the
         755             :     ! intermediate thermodynamic levels, d(var_zm)/dz is multiplied by
         756             :     ! ( K_zt + nu ).  Then, the derivative of the whole mathematical expression
         757             :     ! is taken over the central momentum level, which yields the desired result.
         758             :     !
         759             :     ! ==var_zm================================================= m(k+1)
         760             :     !
         761             :     ! ----------d(var_zm)/dz--(K_zt+nu)------------------------ t(k+1)
         762             :     !
         763             :     ! ==var_zm===================d[(K_zt+nu)*d(var_zm)/dz]/dz== m(k)
         764             :     !
         765             :     ! ----------d(var_zm)/dz--(K_zt+nu)------------------------ t(k)
         766             :     !
         767             :     ! ==var_zm================================================= m(k-1)
         768             :     !
         769             :     ! The vertical indices m(k+1), t(k+1), m(k), t(k), and m(k-1) correspond
         770             :     ! with altitudes zm(k+1), zt(k+1), zm(k), zt(k), and zm(k-1), respectively.
         771             :     ! The letter "t" is used for thermodynamic levels and the letter "m" is used
         772             :     ! for momentum levels.
         773             :     !
         774             :     ! invrs_dzm(k)   = 1 / ( zt(k+1) - zt(k) )
         775             :     ! invrs_dzt(k+1) = 1 / ( zm(k+1) - zm(k) )
         776             :     ! invrs_dzt(k)   = 1 / ( zm(k) - zm(k-1) )
         777             :     !
         778             :     ! Note:  This function only computes the general implicit form:
         779             :     !        + d [ ( K_zt + nu ) * d( var_zm(t+1) )/dz ] / dz.
         780             :     !        For a Crank-Nicholson scheme, the left-hand side result of this
         781             :     !        function will have to be multiplied by (1/2).  For a
         782             :     !        Crank-Nicholson scheme, the right-hand side (explicit) component
         783             :     !        needs to be computed by multiplying the left-hand side results by
         784             :     !        (1/2), reversing the sign on each left-hand side element, and then
         785             :     !        multiplying each element by the appropriate var_zm(t) value from
         786             :     !        the appropriate vertical level.
         787             :     !
         788             :     !
         789             :     ! Boundary Conditions:
         790             :     !
         791             :     ! 1) Zero-flux boundary conditions.
         792             :     !    This function is set up to use zero-flux boundary conditions at both
         793             :     !    the lower boundary level and the upper boundary level.  The flux, F,
         794             :     !    is the amount of var_zm flowing normal through the boundary per unit
         795             :     !    time per unit surface area.  The derivative of the flux effects the
         796             :     !    time-tendency of var_zm, such that:
         797             :     !
         798             :     !    d(var_zm)/dt = -dF/dz.
         799             :     !
         800             :     !    For the 2nd-order eddy-diffusion term, +d[(K_zt+nu)*d(var_zm)/dz]/dz,
         801             :     !    the flux is:
         802             :     !
         803             :     !    F = -(K_zt+nu)*d(var_zm)/dz.
         804             :     !
         805             :     !    In order to have zero-flux boundary conditions, the derivative of
         806             :     !    var_zm, d(var_zm)/dz, needs to equal 0 at both the lower boundary and
         807             :     !    the upper boundary.
         808             :     !
         809             :     !    In order to discretize the lower boundary condition, consider a new
         810             :     !    level outside the model (momentum level 0) just below the lower
         811             :     !    boundary level (momentum level 1).  The value of var_zm at the level
         812             :     !    just outside the model is defined to be the same as the value of var_zm
         813             :     !    at the lower boundary level.  Therefore, the value of d(var_zm)/dz
         814             :     !    between the level just outside the model and the lower boundary level
         815             :     !    is 0, satisfying the zero-flux boundary condition.  The other value for
         816             :     !    d(var_zm)/dz (between momentum level 2 and momentum level 1) is taken
         817             :     !    over the intermediate thermodynamic level (thermodynamic level 2),
         818             :     !    where it is multiplied by the factor ( K_zt(2) + nu ).  Then, the
         819             :     !    derivative of the whole expression is taken over the central momentum
         820             :     !    level.
         821             :     !
         822             :     !    =var_zm(2)============================================ m(2)
         823             :     !
         824             :     !    ----------d(var_zm)/dz==(K_zt(2)+nu)------------------ t(2)
         825             :     !
         826             :     !    =var_zm(1)===============d[(K_zt+nu)*d(var_zm)/dz]/dz= m(1) Boundary
         827             :     !
         828             :     !    ----------[d(var_zm)/dz = 0]-------------------------- t(1)
         829             :     !
         830             :     !    =[var_zm(0) = var_zm(1)]=====(level outside model)==== m(0)
         831             :     !
         832             :     !    The result is dependent only on values of K_zt found at thermodynamic
         833             :     !    level 2 and values of var_zm found at momentum levels 1 and 2.  Thus,
         834             :     !    it only affects 2 diagonals on the left-hand side matrix.
         835             :     !
         836             :     !    The same method can be used to discretize the upper boundary by
         837             :     !    considering a new level outside the model just above the upper boundary
         838             :     !    level.
         839             :     !
         840             :     ! 2) Fixed-point boundary conditions.
         841             :     !    Many equations in the model use fixed-point boundary conditions rather
         842             :     !    than zero-flux boundary conditions.  This means that the value of
         843             :     !    var_zm stays the same over the course of the timestep at the lower
         844             :     !    boundary, as well as at the upper boundary.
         845             :     !
         846             :     !    In order to discretize the boundary conditions for equations requiring
         847             :     !    fixed-point boundary conditions, either:
         848             :     !    a) in the parent subroutine or function (that calls this function),
         849             :     !       loop over all vertical levels from the second-lowest to the
         850             :     !       second-highest, ignoring the boundary levels.  Then set the values
         851             :     !       at the boundary levels in the parent subroutine; or
         852             :     !    b) in the parent subroutine or function, loop over all vertical levels
         853             :     !       and then overwrite the results at the boundary levels.
         854             :     !
         855             :     !    Either way, at the boundary levels, an array with a value of 1 at the
         856             :     !    main diagonal on the left-hand side and with values of 0 at all other
         857             :     !    diagonals on the left-hand side will preserve the right-hand side value
         858             :     !    at that level, thus satisfying the fixed-point boundary conditions.
         859             :     !
         860             :     !
         861             :     ! Conservation Properties:
         862             :     !
         863             :     ! When zero-flux boundary conditions are used, this technique of
         864             :     ! discretizing the eddy diffusion term leads to conservative differencing.
         865             :     ! When conservative differencing is in place, the column totals for each
         866             :     ! column in the left-hand side matrix (for the eddy diffusion term) should
         867             :     ! be equal to 0.  This ensures that the total amount of the quantity var_zm
         868             :     ! over the entire vertical domain is being conserved, meaning that nothing
         869             :     ! is lost due to diffusional effects.
         870             :     !
         871             :     ! To see that this conservation law is satisfied, compute the eddy diffusion
         872             :     ! of var_zm and integrate vertically.  In discretized matrix notation (where
         873             :     ! "i" stands for the matrix column and "j" stands for the matrix row):
         874             :     !
         875             :     !  0 = Sum_j Sum_i ( 1/invrs_dzm )_i
         876             :     !                     ( invrs_dzm * ((K_zt+nu)*invrs_dzt) )_ij (var_zm)_j.
         877             :     !
         878             :     ! The left-hand side matrix, ( invrs_dzm * ((K_zt+nu)*invrs_dzt) )_ij, is
         879             :     ! partially written below.  The sum over i in the above equation removes
         880             :     ! invrs_dzm everywhere from the matrix below.  The sum over j leaves the
         881             :     ! column totals that are desired.
         882             :     !
         883             :     ! Left-hand side matrix contributions from eddy diffusion term; first four
         884             :     ! vertical levels:
         885             :     !
         886             :     !     ---------------------------------------------------------------------->
         887             :     !k=1 | +invrs_dzm(k)          -invrs_dzm(k)                      0
         888             :     !    |   *(K_zt(k+1)+nu)        *(K_zt(k+1)+nu)
         889             :     !    |     *invrs_dzt(k+1)        *invrs_dzt(k+1)
         890             :     !    |
         891             :     !k=2 | -invrs_dzm(k)          +invrs_dzm(k)            -invrs_dzm(k)
         892             :     !    |   *(K_zt(k)+nu)          *[ (K_zt(k+1)+nu)        *(K_zt(k+1)+nu)
         893             :     !    |     *invrs_dzt(k)            *invrs_dzt(k+1)        *invrs_dzt(k+1)
         894             :     !    |                            +(K_zt(k)+nu)
         895             :     !    |                              *invrs_dzt(k) ]
         896             :     !    |
         897             :     !k=3 |          0             -invrs_dzm(k)            +invrs_dzm(k)
         898             :     !    |                          *(K_zt(k)+nu)            *[ (K_zt(k+1)+nu)
         899             :     !    |                            *invrs_dzt(k)              *invrs_dzt(k+1)
         900             :     !    |                                                     +(K_zt(k)+nu)
         901             :     !    |                                                       *invrs_dzt(k) ]
         902             :     !    |
         903             :     !k=4 |          0                       0              -invrs_dzm(k)
         904             :     !    |                                                   *(K_zt(k)+nu)
         905             :     !    |                                                     *invrs_dzt(k)
         906             :     !   \ /
         907             :     !
         908             :     ! Note:  The superdiagonal term from level 3 and both the main diagonal and
         909             :     !        superdiagonal terms from level 4 are not shown on this diagram.
         910             :     !
         911             :     ! Note:  The matrix shown is a tridiagonal matrix.  For a band diagonal
         912             :     !        matrix (with 5 diagonals), there would be an extra row between each
         913             :     !        of the rows shown and an extra column between each of the columns
         914             :     !        shown.  However, for the purposes of the var_zm eddy diffusion
         915             :     !        term, those extra row and column values are all 0, and the
         916             :     !        conservation properties of the matrix aren't effected.
         917             :     !
         918             :     ! If fixed-point boundary conditions are used, the matrix entries at
         919             :     ! level 1 (k=1) read:  1   0   0; which means that conservative differencing
         920             :     ! is not in play.  The total amount of var_zm over the entire vertical
         921             :     ! domain is not being conserved, as amounts of var_zm may be fluxed out
         922             :     ! through the upper boundary or lower boundary through the effects of
         923             :     ! diffusion.
         924             :     !
         925             :     ! Brian Griffin.  April 26, 2008.
         926             :
         927             :     ! References:
         928             :     ! None
         929             :     !-----------------------------------------------------------------------
         930             :
         931             :     use grid_class, only: &
         932             :         grid, & ! Type
         933             :         ddzt    ! Procedure
         934             :
         935             :     use constants_clubb, only: &
         936             :         zero    ! Constant(s)
         937             :
         938             :     use clubb_precision, only: &
         939             :         core_rknd ! Variable(s)
         940             :
         941             :     use model_flags, only: &
         942             :         l_upwind_Kh_dp_term
         943             :
         944             :     implicit none
         945             :
         946             :     ! Constant parameters
         947             :     integer, parameter :: &
         948             :       kp1_mdiag = 1, & ! Momentum superdiagonal index.
         949             :       k_mdiag   = 2, & ! Momentum main diagonal index.
         950             :       km1_mdiag = 3    ! Momentum subdiagonal index.
         951             :
         952             :     ! Input Variables
         953             :     integer, intent(in) :: &
         954             :       nz, &
         955             :       ngrdcol
         956             :
         957             :     type (grid), target, intent(in) :: gr
         958             :
         959             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         960             :       K_zm,           &
         961             :       K_zt,           & ! Coef. of eddy diffusivity at thermo. levels   [m^2/s]
         962             :       rho_ds_zt,      &
         963             :       invrs_rho_ds_zm
         964             :
         965             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
         966             :       nu                ! Background constant coef. of eddy diffusivity [m^2/s]
         967             :
         968             :     ! Return Variable
         969             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(out) :: &
         970             :       lhs     ! LHS coefficient of diffusion term    [1/s]
         971             :
         972             :     ! Local Variable
         973             :     integer :: i, k    ! Vertical level index
         974             :
         975             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz) :: &
         976    71480448 :       lhs_upwind    ! LHS coefficient diffusion term due to upwind  [1/s]
         977             :
         978             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         979    71480448 :       drhoKdz_zm, &
         980    71480448 :       rho_K_zt_nu, &
         981    71480448 :       ddzt_rho_K_zt_nu
         982             :
         983             :     !------------Begin code------------------------------
         984             :
         985             :     !$acc data copyin( gr, gr%invrs_dzt, gr%invrs_dzm, &
         986             :     !$acc              K_zm, K_zt, rho_ds_zt, invrs_rho_ds_zm, nu ) &
         987             :     !$acc     copyout( lhs ) &
         988             :     !$acc      create( lhs_upwind, drhoKdz_zm, rho_K_zt_nu, ddzt_rho_K_zt_nu )
         989             :
         990             :     !$acc parallel loop gang vector collapse(2) default(present)
         991  3073659264 :     do k = 1, nz
         992 50761923264 :       do i = 1, ngrdcol
         993 50726183040 :         rho_K_zt_nu(i,k) = rho_ds_zt(i,k) * ( K_zt(i,k) + nu(i) )
         994             :       end do
         995             :     end do
         996             :     !$acc end parallel loop
         997             :
         998             :     ! calculate the dKh_zm/dz
         999    35740224 :     ddzt_rho_K_zt_nu = ddzt( nz, ngrdcol, gr, rho_K_zt_nu )
        1000             :
        1001             :     !$acc parallel loop gang vector collapse(2) default(present)
        1002  3073659264 :     do k = 1, nz
        1003 50761923264 :       do i = 1, ngrdcol
        1004 50726183040 :         drhoKdz_zm(i,k) = - invrs_rho_ds_zm(i,k) * ddzt_rho_K_zt_nu(i,k)
        1005             :       end do
        1006             :     end do
        1007             :     !$acc end parallel loop
        1008             :
        1009             :     ! extra terms with upwind scheme
        1010             :     ! k = 1 (bottom level); lowere boundary level
        1011             :     !$acc parallel loop gang vector default(present)
        1012   596778624 :     do i = 1, ngrdcol
        1013   561038400 :       lhs_upwind(kp1_mdiag,i,1) = + min( drhoKdz_zm(i,1) , zero ) * gr%invrs_dzt(i,2)
        1014   561038400 :       lhs_upwind(k_mdiag,i,1)   = - min( drhoKdz_zm(i,1) , zero ) * gr%invrs_dzt(i,2)
        1015   596778624 :       lhs_upwind(km1_mdiag,i,1) = zero
        1016             :     end do
        1017             :     !$acc end parallel loop
        1018             :
        1019             :     ! Most of the interior model; normal conditions.
        1020             :     !$acc parallel loop gang vector collapse(2) default(present)
        1021  3002178816 :     do k = 2, nz-1, 1
        1022 49568366016 :       do i = 1, ngrdcol
        1023             :         ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
        1024 46566187200 :         lhs_upwind(kp1_mdiag,i,k) = + min( drhoKdz_zm(i,k) , zero ) * gr%invrs_dzt(i,k+1)
        1025             :         ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
        1026           0 :         lhs_upwind(k_mdiag,i,k)   = - min( drhoKdz_zm(i,k) , zero ) * gr%invrs_dzt(i,k+1) &
        1027 46566187200 :                                     + max( drhoKdz_zm(i,k) , zero ) * gr%invrs_dzt(i,k)
        1028             :         ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
        1029 49532625792 :         lhs_upwind(km1_mdiag,i,k) = - max( drhoKdz_zm(i,k) , zero ) * gr%invrs_dzt(i,k)
        1030             :       end do
        1031             :     end do
        1032             :     !$acc end parallel loop
        1033             :
        1034             :     ! k = nz (top level); upper boundary level.
        1035             :     ! Only relevant if zero-flux boundary conditions are used.
        1036             :     !$acc parallel loop gang vector default(present)
        1037   596778624 :     do i = 1, ngrdcol
        1038   561038400 :       lhs_upwind(kp1_mdiag,i,nz) = zero
        1039   561038400 :       lhs_upwind(k_mdiag,i,nz)   = + max( drhoKdz_zm(i,nz) , zero ) * gr%invrs_dzt(i,nz)
        1040   596778624 :       lhs_upwind(km1_mdiag,i,nz) = - max( drhoKdz_zm(i,nz) , zero ) * gr%invrs_dzt(i,nz)
        1041             :     end do
        1042             :     !$acc end parallel loop
        1043             :
        1044             :     ! k = 1; lower boundary level at surface.
        1045             :     ! Only relevant if zero-flux boundary conditions are used.
        1046             :     ! These k=1 lines currently do not have any effect on model results.
        1047             :     ! These k=1 level of this "lhs" array is not fed into
        1048             :     ! the final LHS matrix that will be used to solve for the next timestep.
        1049             :
        1050             :     if ( .not. l_upwind_Kh_dp_term ) then
        1051             :
        1052             :       !$acc parallel loop gang vector default(present)
        1053   596778624 :       do i = 1, ngrdcol
        1054             :         ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
        1055  1122076800 :         lhs(kp1_mdiag,i,1) = - gr%invrs_dzm(i,1) * invrs_rho_ds_zm(i,1) &
        1056  1122076800 :                                 * ( K_zt(i,2) + nu(i) ) * rho_ds_zt(i,2) * gr%invrs_dzt(i,2)
        1057             :
        1058             :         ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
        1059           0 :         lhs(k_mdiag,i,1)   = + gr%invrs_dzm(i,1) * invrs_rho_ds_zm(i,1) &
        1060   561038400 :                                   * ( K_zt(i,2) + nu(i) ) * rho_ds_zt(i,2) * gr%invrs_dzt(i,2)
        1061             :
        1062             :         ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
        1063   596778624 :         lhs(km1_mdiag,i,1) = zero
        1064             :       end do
        1065             :       !$acc end parallel loop
        1066             :
        1067             :     else
        1068             :
        1069             :       !$acc parallel loop gang vector default(present)
        1070             :       do i = 1, ngrdcol
        1071             :         ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
        1072             :         lhs(kp1_mdiag,i,1) = - gr%invrs_dzm(i,1) * ( K_zm(i,1) + nu(i) ) * gr%invrs_dzt(i,2) &
        1073             :                            + lhs_upwind(kp1_mdiag,i,1)
        1074             :
        1075             :         ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
        1076             :         lhs(k_mdiag,i,1)   = + gr%invrs_dzm(i,1) * ( K_zm(i,1) + nu(i) ) * gr%invrs_dzt(i,2) &
        1077             :                            + lhs_upwind(k_mdiag,i,1)
        1078             :
        1079             :         ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
        1080             :         lhs(km1_mdiag,i,1) =   zero + lhs_upwind(km1_mdiag,i,1)
        1081             :       end do
        1082             :       !$acc end parallel loop
        1083             :
        1084             :     end if
        1085             :
        1086             :     ! Most of the interior model; normal conditions.
        1087             :     if ( .not. l_upwind_Kh_dp_term ) then
        1088             :
        1089             :       !$acc parallel loop gang vector collapse(2) default(present)
        1090  3002178816 :       do k = 2, nz-1, 1
        1091 49568366016 :         do i = 1, ngrdcol
        1092             :           ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
        1093 93132374400 :           lhs(kp1_mdiag,i,k) &
        1094           0 :           = - gr%invrs_dzm(i,k) * invrs_rho_ds_zm(i,k) &
        1095 >13969*10^7 :                            * ( K_zt(i,k+1) + nu(i) ) * rho_ds_zt(i,k+1) * gr%invrs_dzt(i,k+1)
        1096             :
        1097             :           ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
        1098             :           lhs(k_mdiag,i,k) &
        1099           0 :           = + gr%invrs_dzm(i,k) * invrs_rho_ds_zm(i,k) &
        1100           0 :                            * ( ( K_zt(i,k+1) + nu(i) ) * rho_ds_zt(i,k+1) * gr%invrs_dzt(i,k+1)  &
        1101 46566187200 :                                + ( K_zt(i,k) + nu(i) ) * rho_ds_zt(i,k) * gr%invrs_dzt(i,k) )
        1102             :
        1103             :           ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
        1104             :           lhs(km1_mdiag,i,k) &
        1105           0 :           = - gr%invrs_dzm(i,k) * invrs_rho_ds_zm(i,k) &
        1106 49532625792 :                            * ( K_zt(i,k) + nu(i) ) * rho_ds_zt(i,k) * gr%invrs_dzt(i,k)
        1107             :         end do
        1108             :       end do
        1109             :       !$acc end parallel loop
        1110             :
        1111             :     else
        1112             :
        1113             :       !$acc parallel loop gang vector collapse(2) default(present)
        1114             :       do k = 2, nz-1, 1
        1115             :         do i = 1, ngrdcol
        1116             :
        1117             :           ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
        1118             :           lhs(kp1_mdiag,i,k) &
        1119             :           = - gr%invrs_dzm(i,k) * ( K_zm(i,k) + nu(i) ) * gr%invrs_dzt(i,k+1) &
        1120             :             + lhs_upwind(kp1_mdiag,i,k)
        1121             :
        1122             :           ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
        1123             :           lhs(k_mdiag,i,k) &
        1124             :           = + gr%invrs_dzm(i,k) * ( K_zm(i,k) + nu(i) ) * ( gr%invrs_dzt(i,k+1) + gr%invrs_dzt(i,k) ) &
        1125             :             + lhs_upwind(k_mdiag,i,k)
        1126             :
        1127             :           ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
        1128             :           lhs(km1_mdiag,i,k) &
        1129             :           = - gr%invrs_dzm(i,k) * ( K_zm(i,k) + nu(i) ) * gr%invrs_dzt(i,k) &
        1130             :             + lhs_upwind(km1_mdiag,i,k)
        1131             :
        1132             :         end do
        1133             :       end do
        1134             :       !$acc end parallel loop
        1135             :
        1136             :     end if
        1137             :
        1138             :     ! k = nz (top level); upper boundary level.
        1139             :     ! Only relevant if zero-flux boundary conditions are used.
        1140             :
        1141             :     if ( .not. l_upwind_Kh_dp_term ) then
        1142             :
        1143             :       !$acc parallel loop gang vector default(present)
        1144   596778624 :       do i = 1, ngrdcol
        1145             :         ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
        1146   561038400 :         lhs(kp1_mdiag,i,nz) = zero
        1147             :
        1148             :         ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
        1149             :         lhs(k_mdiag,i,nz) &
        1150           0 :         = + gr%invrs_dzm(i,nz) * invrs_rho_ds_zm(i,nz) &
        1151   561038400 :              * ( K_zt(i,nz) + nu(i) ) * rho_ds_zt(i,nz) * gr%invrs_dzt(i,nz)
        1152             :
        1153             :         ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
        1154             :         lhs(km1_mdiag,i,nz) &
        1155           0 :         = - gr%invrs_dzm(i,nz) * invrs_rho_ds_zm(i,nz) &
        1156   596778624 :              * ( K_zt(i,nz) + nu(i) ) * rho_ds_zt(i,nz) * gr%invrs_dzt(i,nz)
        1157             :       end do
        1158             :       !$acc end parallel loop
        1159             :
        1160             :     else
        1161             :
        1162             :       !$acc parallel loop gang vector default(present)
        1163             :       do i = 1, ngrdcol
        1164             :         ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
        1165             :         lhs(kp1_mdiag,i,nz) =  zero &
        1166             :                               + lhs_upwind(kp1_mdiag,i,nz)
        1167             :
        1168             :         ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
        1169             :         lhs(k_mdiag,i,nz) &
        1170             :         = + gr%invrs_dzm(i,nz) * ( K_zm(i,nz) + nu(i) ) * gr%invrs_dzt(i,nz) &
        1171             :           + lhs_upwind(k_mdiag,i,nz)
        1172             :
        1173             :         ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
        1174             :         lhs(km1_mdiag,i,nz) &
        1175             :         = - gr%invrs_dzm(i,nz) * ( K_zm(i,nz) + nu(i) ) * gr%invrs_dzt(i,nz) &
        1176             :           + lhs_upwind(km1_mdiag,i,nz)
        1177             :       end do
        1178             :       !$acc end parallel loop
        1179             :
        1180             :     end if
        1181             :
        1182             :     !$acc end data
        1183             :
        1184    35740224 :     return
        1185             :
        1186             :   end subroutine diffusion_zm_lhs
        1187             :
        1188             :   !=============================================================================
        1189             :
        1190             : end module diffusion

Generated by: LCOV version 1.14