LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - advance_windm_edsclrm_module.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 143 | 489 | 29.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 4 | 6 | 66.7 %

              Line data    Source code
           1             :
           2             : module advance_windm_edsclrm_module
           3             :
           4             :   implicit none
           5             :
           6             :   private ! Set Default Scope
           7             :
           8             :   public :: advance_windm_edsclrm
           9             :
          10             :   private :: windm_edsclrm_solve, &
          11             :     compute_uv_tndcy,  &
          12             :     windm_edsclrm_lhs, &
          13             :     windm_edsclrm_rhs
          14             :
          15             :
          16             :   ! Private named constants to avoid string comparisons
          17             :   integer, parameter, private :: &
          18             :     windm_edsclrm_um = 1, &     ! Named constant to handle um solves
          19             :     windm_edsclrm_vm = 2, &     ! Named constant to handle vm solves
          20             :     windm_edsclrm_scalar = 3, & ! Named constant to handle scalar solves
          21             :     clip_upwp = 10, &           ! Named constant for upwp clipping
          22             :                                 ! NOTE: This must be the same as the clip_upwp
          23             :                                 ! declared in clip_explicit!
          24             :     clip_vpwp = 11              ! Named constant for vpwp clipping
          25             :                                 ! NOTE: This must be the same as the clip_vpwp
          26             :                                 ! declared in clip_explicit!
          27             :
          28             :   contains
          29             :
          30             :   !=============================================================================
          31     8935056 :   subroutine advance_windm_edsclrm( nz, ngrdcol, edsclr_dim, gr, dt, &
          32     8935056 :                                     wm_zt, Km_zm, Kmh_zm, &
          33     8935056 :                                     ug, vg, um_ref, vm_ref, &
          34     8935056 :                                     wp2, up2, vp2, um_forcing, vm_forcing, &
          35     8935056 :                                     edsclrm_forcing, &
          36     8935056 :                                     rho_ds_zm, invrs_rho_ds_zt, &
          37     8935056 :                                     fcor, l_implemented, &
          38             :                                     nu_vert_res_dep, &
          39             :                                     tridiag_solve_method, &
          40             :                                     l_predict_upwp_vpwp, &
          41             :                                     l_upwind_xm_ma, &
          42             :                                     l_uv_nudge, &
          43             :                                     l_tke_aniso, &
          44             :                                     l_lmm_stepping, &
          45             :                                     l_linearize_pbl_winds, &
          46             :                                     order_xp2_xpyp, order_wp2_wp3, order_windm, &
          47             :                                     stats_metadata, &
          48     8935056 :                                     stats_zt, stats_zm, stats_sfc, &
          49     8935056 :                                     um, vm, edsclrm, &
          50     8935056 :                                     upwp, vpwp, wpedsclrp, &
          51     8935056 :                                     um_pert, vm_pert, upwp_pert, vpwp_pert )
          52             :
          53             :     ! Description:
          54             :     ! Solves for both mean horizontal wind components, um and vm, and for the
          55             :     ! eddy-scalars (passive scalars that don't use the high-order closure).
          56             :
          57             :     ! Uses the LAPACK tridiagonal solver subroutine with 2 + # of scalar(s)
          58             :     ! back substitutions (since the left hand side matrix is the same for all
          59             :     ! input variables).
          60             :
          61             :     ! References:
          62             :     ! Eqn. 8 & 9 on p. 3545 of
          63             :     ! ``A PDF-Based Model for Boundary Layer Clouds. Part I:
          64             :     ! Method and Model Description'' Golaz, et al. (2002)
          65             :     ! JAS, Vol. 59, pp. 3540--3551.
          66             :     !-----------------------------------------------------------------------
          67             :
          68             :     use grid_class, only:  &
          69             :         grid, &  ! Type
          70             :         zm2zt
          71             :
          72             :     use parameters_model, only:  &
          73             :         ts_nudge ! Variable(s)
          74             :
          75             :     use parameters_tunable, only: &
          76             :         nu_vertical_res_dep    ! Type(s)
          77             :
          78             :     use clubb_precision, only:  &
          79             :         core_rknd ! Variable(s)
          80             :
          81             :     use stats_type_utilities, only: &
          82             :         stat_begin_update, & ! Subroutines
          83             :         stat_end_update, &
          84             :         stat_update_var
          85             :
          86             :     use stats_variables, only: &
          87             :         stats_metadata_type
          88             :
          89             :     use clip_explicit, only:  &
          90             :         clip_covar  ! Procedure(s)
          91             :
          92             :     use error_code, only: &
          93             :         clubb_at_least_debug_level,  & ! Procedure
          94             :         err_code,                    & ! Error Indicator
          95             :         clubb_fatal_error              ! Constant
          96             :
          97             :     use constants_clubb, only:  &
          98             :         one_half, & ! Constant(s)
          99             :         zero,     &
         100             :         fstderr,  &
         101             :         eps
         102             :
         103             :     use sponge_layer_damping, only: &
         104             :         uv_sponge_damp_settings, &
         105             :         uv_sponge_damp_profile, &
         106             :         sponge_damp_xm     ! Procedure(s)
         107             :
         108             :     use stats_type, only: stats ! Type
         109             :
         110             :     use diffusion, only:  &
         111             :         diffusion_zt_lhs   ! Procedure(s)
         112             :
         113             :     use mean_adv, only: &
         114             :         term_ma_zt_lhs    ! Procedures
         115             :
         116             :     use model_flags, only: &
         117             :         l_upwind_Kh_dp_term
         118             :
         119             :     use advance_helper_module, only: &
         120             :         calc_xpwp
         121             :
         122             :     implicit none
         123             :
         124             :     ! ------------------------ Input Variables ------------------------
         125             :     integer, intent(in) :: &
         126             :       nz, &
         127             :       ngrdcol, &
         128             :       edsclr_dim
         129             :
         130             :     type (grid), target, intent(in) :: &
         131             :       gr
         132             :
         133             :     real( kind = core_rknd ), intent(in) ::  &
         134             :       dt                 ! Model timestep                             [s]
         135             :
         136             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         137             :       wm_zt,           & ! w wind component on thermodynamic levels    [m/s]
         138             :       Km_zm,           & ! Eddy diffusivity of winds on momentum levs. [m^2/s]
         139             :       Kmh_zm,          & ! Eddy diffusivity of themo on momentum levs. [m^s/s]
         140             :       ug,              & ! u (west-to-east) geostrophic wind comp.     [m/s]
         141             :       vg,              & ! v (south-to-north) geostrophic wind comp.   [m/s]
         142             :       um_ref,          & ! Reference u wind component for nudging      [m/s]
         143             :       vm_ref,          & ! Reference v wind component for nudging      [m/s]
         144             :       wp2,             & ! w'^2 (momentum levels)                      [m^2/s^2]
         145             :       up2,             & ! u'^2 (momentum levels)                      [m^2/s^2]
         146             :       vp2,             & ! v'^2 (momentum levels)                      [m^2/s^2]
         147             :       um_forcing,      & ! u forcing                                   [m/s/s]
         148             :       vm_forcing,      & ! v forcing                                   [m/s/s]
         149             :       rho_ds_zm,       & ! Dry, static density on momentum levels      [kg/m^3]
         150             :       invrs_rho_ds_zt    ! Inv. dry, static density at thermo. levels  [m^3/kg]
         151             :
         152             :     real( kind = core_rknd ), dimension(ngrdcol,nz,edsclr_dim), intent(in) ::  &
         153             :       edsclrm_forcing  ! Eddy scalar large-scale forcing        [{units vary}/s]
         154             :
         155             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
         156             :       fcor           ! Coriolis parameter                            [s^-1]
         157             :
         158             :     logical, intent(in) ::  &
         159             :       l_implemented  ! Flag for CLUBB being implemented in a larger model.
         160             :
         161             :     type(nu_vertical_res_dep), intent(in) :: &
         162             :       nu_vert_res_dep    ! Vertical resolution dependent nu values
         163             :
         164             :     integer, intent(in) :: &
         165             :       tridiag_solve_method  ! Specifier for method to solve tridiagonal systems
         166             :
         167             :     logical, intent(in) :: &
         168             :       l_predict_upwp_vpwp,   & ! Flag to predict <u'w'> and <v'w'> along with <u> and <v> alongside
         169             :                                ! the advancement of <rt>, <w'rt'>, <thl>, <wpthlp>, <sclr>, and
         170             :                                ! <w'sclr'> in subroutine advance_xm_wpxp.  Otherwise, <u'w'> and
         171             :                                ! <v'w'> are still approximated by eddy diffusivity when <u> and <v>
         172             :                                ! are advanced in subroutine advance_windm_edsclrm.
         173             :       l_upwind_xm_ma,        & ! This flag determines whether we want to use an upwind differencing
         174             :                                ! approximation rather than a centered differencing for turbulent or
         175             :                                ! mean advection terms. It affects rtm, thlm, sclrm, um and vm.
         176             :       l_uv_nudge,            & ! For wind speed nudging
         177             :       l_tke_aniso,           & ! For anisotropic turbulent kinetic energy, i.e. TKE = 1/2
         178             :                                ! (u'^2 + v'^2 + w'^2)
         179             :       l_lmm_stepping,        & ! Apply Linear Multistep Method (LMM) Stepping
         180             :       l_linearize_pbl_winds    ! Flag (used by E3SM) to linearize PBL winds
         181             :
         182             :     integer, intent(in) :: &
         183             :       order_xp2_xpyp, &
         184             :       order_wp2_wp3, &
         185             :       order_windm
         186             :
         187             :     type (stats_metadata_type), intent(in) :: &
         188             :       stats_metadata
         189             :
         190             :     ! ------------------------ Input/Output Variables ------------------------
         191             :     type (stats), dimension(ngrdcol), target, intent(inout) :: &
         192             :       stats_zt, &
         193             :       stats_zm, &
         194             :       stats_sfc
         195             :
         196             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) ::  &
         197             :       um,   & ! Mean u (west-to-east) wind component         [m/s]
         198             :       vm,   & ! Mean v (south-to-north) wind component       [m/s]
         199             :       upwp, & ! <u'w'> (momentum levels)                     [m^2/s^2]
         200             :       vpwp    ! <v'w'> (momentum levels)                     [m^2/s^2]
         201             :
         202             :     ! Input/Output Variable for eddy-scalars
         203             :     real( kind = core_rknd ), dimension(ngrdcol,nz,edsclr_dim), intent(inout) ::  &
         204             :       edsclrm        ! Mean eddy scalar quantity             [units vary]
         205             :
         206             :     ! Output Variable for eddy-scalars
         207             :     real( kind = core_rknd ), dimension(ngrdcol,nz,edsclr_dim), intent(inout) ::  &
         208             :       wpedsclrp      ! w'edsclr' (momentum levels)           [m/s {units vary}]
         209             :
         210             :     ! Variables used to track perturbed version of winds.
         211             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
         212             :       um_pert,   & ! perturbed <u>    [m/s]
         213             :       vm_pert,   & ! perturbed <v>    [m/s]
         214             :       upwp_pert, & ! perturbed <u'w'> [m^2/s^2]
         215             :       vpwp_pert    ! perturbed <v'w'> [m^2/s^2]
         216             :
         217             :     ! ------------------------ Local Variables ------------------------
         218             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         219    17870112 :       um_old, & ! Saved value of mean u (west-to-east) wind component    [m/s]
         220    17870112 :       vm_old    ! Saved value of Mean v (south-to-north) wind component  [m/s]
         221             :
         222             :     real( kind = core_rknd ), dimension(ngrdcol,nz,edsclr_dim) ::  &
         223    17870112 :       edsclrm_old    ! Saved value of mean eddy scalar quantity   [units vary]
         224             :
         225             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         226    17870112 :       um_tndcy,    & ! u wind component tendency                    [m/s^2]
         227    17870112 :       vm_tndcy       ! v wind component tendency                    [m/s^2]
         228             :
         229             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         230    17870112 :       upwp_chnge,  & ! Net change of u'w' due to clipping           [m^2/s^2]
         231    17870112 :       vpwp_chnge     ! Net change of v'w' due to clipping           [m^2/s^2]
         232             :
         233             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz) :: &
         234    17870112 :       lhs ! The implicit part of the tridiagonal matrix             [units vary]
         235             :
         236             :     real( kind = core_rknd ), dimension(ngrdcol,nz,max(2,edsclr_dim)) :: &
         237    17870112 :       rhs,     &! The explicit part of the tridiagonal matrix       [units vary]
         238    17870112 :       solution  ! The solution to the tridiagonal matrix            [units vary]
         239             :
         240             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         241    17870112 :       wind_speed,      & ! wind speed; sqrt(u^2 + v^2)              [m/s]
         242    17870112 :       wind_speed_pert    ! perturbed wind speed; sprt(u^2 + v^2)    [m/s]
         243             :
         244             :     real( kind = core_rknd ), dimension(ngrdcol) :: &
         245    17870112 :       u_star_sqd,      & ! Surface friction velocity, u_star, squared      [m/s]
         246    17870112 :       u_star_sqd_pert    ! perturbed u_star, squared                       [m/s]
         247             :
         248             :     logical :: &
         249             :       l_imp_sfc_momentum_flux  ! Flag for implicit momentum surface fluxes.
         250             :
         251             :     integer :: nrhs  ! Number of right hand side terms
         252             :
         253             :     integer :: i, k, edsclr, j  ! Array index
         254             :
         255             :     logical :: l_first_clip_ts, l_last_clip_ts ! flags for clip_covar
         256             :
         257             :     real( kind = core_rknd ), dimension(ngrdcol) :: &
         258    17870112 :       nu_zero
         259             :
         260             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz) :: &
         261    17870112 :       lhs_diff, & ! LHS diffustion terms
         262    17870112 :       lhs_ma_zt   ! LHS mean advection terms
         263             :
         264             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         265    17870112 :       Km_zt,          & ! Eddy diffusivity of winds on momentum levs. [m^2/s]
         266    17870112 :       Kmh_zt,         & ! Eddy diffusivity of themo on momentum levs. [m^s/s]
         267    17870112 :       Km_zm_p_nu10,   & ! Km_zm plus nu_vert_res_dep%nu10
         268    17870112 :       xpwp              ! x'w' for arbitrary x
         269             :
         270             :     integer, parameter :: &
         271             :       kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
         272             :       k_tdiag   = 2, & ! Thermodynamic main diagonal index.
         273             :       km1_tdiag = 3    ! Thermodynamic subdiagonal index.
         274             :
         275             :     ! Whether perturbed winds are being solved.
         276             :     logical :: l_perturbed_wind
         277             :
         278             :     real( kind = core_rknd ), dimension(nz) ::  &
         279    17870112 :       tmp_in
         280             :
         281             :     ! ------------------------ Begin Code ------------------------
         282             :
         283             :     !$acc enter data create( um_old, vm_old, um_tndcy, vm_tndcy, &
         284             :     !$acc                    upwp_chnge, vpwp_chnge, lhs, rhs, solution, wind_speed, &
         285             :     !$acc                    wind_speed_pert, u_star_sqd, u_star_sqd_pert, &
         286             :     !$acc                    nu_zero, lhs_diff, lhs_ma_zt, Km_zt, Kmh_zt, &
         287             :     !$acc                    Km_zm_p_nu10, xpwp )
         288             :
         289             :     !$acc enter data if( edsclr_dim > 0) create( edsclrm_old )
         290             :
         291             :     !$acc parallel loop gang vector default(present)
         292   149194656 :     do i = 1, ngrdcol
         293   149194656 :       nu_zero(i) = zero
         294             :     end do
         295             :     !$acc end parallel loop
         296             :
         297             :     !$acc parallel loop gang vector collapse(2) default(present)
         298   768414816 :     do k = 1, nz
         299 12690480816 :       do i = 1, ngrdcol
         300 12681545760 :         Km_zm_p_nu10(i,k) = Km_zm(i,k) + nu_vert_res_dep%nu10(i)
         301             :       end do
         302             :     end do
         303             :     !$acc end parallel loop
         304             :
         305     8935056 :     l_perturbed_wind = ( .not. l_predict_upwp_vpwp ) .and. l_linearize_pbl_winds
         306             :
         307     8935056 :     if ( .not. l_implemented ) then
         308             :       call term_ma_zt_lhs( nz, ngrdcol, wm_zt, gr%weights_zt2zm, & ! intent(in)
         309             :                            gr%invrs_dzt, gr%invrs_dzm,    & ! intent(in)
         310             :                            l_upwind_xm_ma,                          & ! intent(in)
         311           0 :                            lhs_ma_zt )                         ! intent(out)
         312             :     else
         313             :       !$acc parallel loop gang vector collapse(3) default(present)
         314   768414816 :       do k = 1, nz
         315 12690480816 :         do i = 1, ngrdcol
         316 48447743760 :           do j = 1, 3
         317 47688264000 :             lhs_ma_zt(j,i,k) = zero
         318             :           end do
         319             :         end do
         320             :       end do
         321             :       !$acc end parallel loop
         322             :     end if
         323             :
         324     8935056 :     if ( .not. l_predict_upwp_vpwp ) then
         325             :
         326           0 :       Km_zt(:,:) = zm2zt( nz, ngrdcol, gr, Km_zm(:,:), zero )
         327             :
         328             :       ! Calculate diffusion terms
         329             :       call diffusion_zt_lhs( nz, ngrdcol, gr, Km_zm, Km_zt, nu_vert_res_dep%nu10, & ! In
         330             :                              invrs_rho_ds_zt, rho_ds_zm,                          & ! In
         331           0 :                              lhs_diff )                                             ! Out
         332             :
         333             :       ! The lower boundary condition needs to be applied here at level 2.
         334             :       if ( .not. l_upwind_Kh_dp_term ) then
         335             :         !$acc parallel loop gang vector default(present)
         336           0 :         do i = 1, ngrdcol
         337             :           ! Thermodynamic superdiagonal: [ x xm(k+1,<t+1>) ]
         338           0 :           lhs_diff(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
         339           0 :                                     * ( Km_zm(i,2) + nu_vert_res_dep%nu10(i) ) &
         340           0 :                                       * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
         341             :
         342             :           ! Thermodynamic main diagonal: [ x xm(k,<t+1>) ]
         343           0 :           lhs_diff(k_tdiag,i,2) = + gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
         344             :                                   * ( Km_zm(i,2) + nu_vert_res_dep%nu10(i) ) &
         345           0 :                                     * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
         346             :
         347             :           ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
         348           0 :           lhs_diff(km1_tdiag,i,2) = zero
         349             :         end do
         350             :         !$acc end parallel loop
         351             :       else
         352             :         !$acc parallel loop gang vector default(present)
         353             :         do i = 1, ngrdcol
         354             :           ! Thermodynamic superdiagonal: [ x xm(k+1,<t+1>) ]
         355             :           lhs_diff(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) &
         356             :                                       * ( Km_zt(i,2) + nu_vert_res_dep%nu10(i) ) &
         357             :                                         * gr%invrs_dzm(i,2)
         358             :
         359             :           ! Thermodynamic main diagonal: [ x xm(k,<t+1>) ]
         360             :           lhs_diff(k_tdiag,i,2) = + gr%invrs_dzt(i,2) &
         361             :                                     * ( Km_zt(i,2) + nu_vert_res_dep%nu10(i) ) &
         362             :                                       * gr%invrs_dzm(i,2)
         363             :
         364             :           ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
         365             :           lhs_diff(km1_tdiag,i,2) = zero
         366             :         end do
         367             :         !$acc end parallel loop
         368             :       end if
         369             :
         370           0 :       if ( l_lmm_stepping ) then
         371             :         !$acc parallel loop gang vector collapse(2) default(present)
         372           0 :         do k = 1, nz
         373           0 :           do i = 1, ngrdcol
         374           0 :             um_old(i,k) = um(i,k)
         375           0 :             vm_old(i,k) = vm(i,k)
         376             :           end do
         377             :         end do
         378             :         !$acc end parallel loop
         379             :       end if ! l_lmm_stepping
         380             :
         381             :       !----------------------------------------------------------------
         382             :       ! Prepare tridiagonal system for horizontal winds, um and vm
         383             :       !----------------------------------------------------------------
         384             :
         385             :       ! Compute Coriolis, geostrophic, and other prescribed wind forcings for um.
         386             :       call compute_uv_tndcy( nz, ngrdcol, windm_edsclrm_um, & ! intent(in)
         387             :                              fcor, vm, vg,                  & ! intent(in)
         388             :                              um_forcing, l_implemented,     & ! intent(in)
         389             :                              stats_metadata,                & ! intent(in)
         390             :                              stats_zt,                      & ! intent(inout)
         391           0 :                              um_tndcy )                       ! intent(out)
         392             :
         393             :       ! Compute Coriolis, geostrophic, and other prescribed wind forcings for vm.
         394             :       call compute_uv_tndcy( nz, ngrdcol, windm_edsclrm_vm, & ! intent(in)
         395             :                              fcor, um, ug,                  & ! intent(in)
         396             :                              vm_forcing, l_implemented,     & ! intent(in)
         397             :                              stats_metadata,                & ! intent(in)
         398             :                              stats_zt,                      & ! intent(inout)
         399           0 :                              vm_tndcy )                       ! intent(out)
         400             :
         401             :       ! Momentum surface fluxes, u'w'|_sfc and v'w'|_sfc, are applied through
         402             :       ! an implicit method, such that:
         403             :       !    x'w'|_sfc = - ( u_star(t)^2 / wind_speed(t) ) * xm(t+1).
         404           0 :       l_imp_sfc_momentum_flux = .true.
         405             :
         406             :       ! Compute wind speed (use threshold "eps" to prevent divide-by-zero
         407             :       ! error).
         408             :       !$acc parallel loop gang vector collapse(2) default(present)
         409           0 :       do k = 1, nz
         410           0 :         do i = 1, ngrdcol
         411           0 :           wind_speed(i,k) = max( sqrt( um(i,k)**2 + vm(i,k)**2 ), eps )
         412             :         end do
         413             :       end do
         414             :       !$acc end parallel loop
         415             :
         416             :       ! Compute u_star_sqd according to the definition of u_star.
         417             :       !$acc parallel loop gang vector default(present)
         418           0 :       do i = 1, ngrdcol
         419           0 :         u_star_sqd(i) = sqrt( upwp(i,1)**2 + vpwp(i,1)**2 )
         420             :       end do
         421             :       !$acc end parallel loop
         422             :
         423             :       ! Compute the explicit portion of the um equation.
         424             :       ! Build the right-hand side vector.
         425             :       call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_um, dt,  & ! intent(in)
         426             :                               lhs_diff, um, um_tndcy,                 & ! intent(in)
         427             :                               rho_ds_zm, invrs_rho_ds_zt,             & ! intent(in)
         428             :                               l_imp_sfc_momentum_flux, upwp(:,1),     & ! intent(in)
         429             :                               stats_metadata,                         & ! intent(in)
         430             :                               stats_zt,                               & ! intent(inout)
         431           0 :                               rhs(:,:,windm_edsclrm_um) )               ! intent(out)
         432             :
         433             :       ! Compute the explicit portion of the vm equation.
         434             :       ! Build the right-hand side vector.
         435             :       call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_vm, dt,  & ! intent(in)
         436             :                               lhs_diff, vm, vm_tndcy,                 & ! intent(in)
         437             :                               rho_ds_zm, invrs_rho_ds_zt,             & ! intent(in)
         438             :                               l_imp_sfc_momentum_flux, vpwp(:,1),     & ! intent(in)
         439             :                               stats_metadata,                         & ! intent(in)
         440             :                               stats_zt,                               & ! intent(inout)
         441           0 :                               rhs(:,:,windm_edsclrm_vm) )               ! intent(out)
         442             :
         443             :       ! Store momentum flux (explicit component)
         444             :
         445             :       ! The surface flux, x'w'(1) = x'w'|_sfc, is set elsewhere in the model.
         446             :       ! upwp(1) = upwp_sfc
         447             :       ! vpwp(1) = vpwp_sfc
         448             :
         449             :       call calc_xpwp( nz, ngrdcol, gr, &
         450             :                       Km_zm_p_nu10, um, &
         451           0 :                       xpwp )
         452             :
         453             :       ! Solve for x'w' at all intermediate model levels.
         454             :       ! A Crank-Nicholson timestep is used.
         455             :       !$acc parallel loop gang vector collapse(2) default(present)
         456           0 :       do k = 2, nz-1
         457           0 :         do i = 1, ngrdcol
         458           0 :           upwp(i,k) = -one_half * xpwp(i,k)
         459             :         end do
         460             :       end do
         461             :       !$acc end parallel loop
         462             :
         463             :       call calc_xpwp( nz, ngrdcol, gr, &
         464             :                       Km_zm_p_nu10, vm, &
         465           0 :                       xpwp )
         466             :
         467             :       !$acc parallel loop gang vector collapse(2) default(present)
         468           0 :       do k = 2, nz-1
         469           0 :         do i = 1, ngrdcol
         470           0 :           vpwp(i,k) = -one_half * xpwp(i,k)
         471             :         end do
         472             :       end do
         473             :       !$acc end parallel loop
         474             :
         475             :       ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
         476             :       ! means that x'w' at the top model level is 0,
         477             :       ! since x'w' = - K_zm * d(xm)/dz.
         478             :       !$acc parallel loop gang vector default(present)
         479           0 :       do i = 1, ngrdcol
         480           0 :         upwp(i,nz) = zero
         481           0 :         vpwp(i,nz) = zero
         482             :       end do
         483             :       !$acc end parallel loop
         484             :
         485             :       ! Compute the implicit portion of the um and vm equations.
         486             :       ! Build the left-hand side matrix.
         487             :       call windm_edsclrm_lhs( nz, ngrdcol, gr, dt,                    & ! intent(in)
         488             :                               lhs_ma_zt, lhs_diff,                    & ! intent(in)
         489             :                               wind_speed, u_star_sqd,                 & ! intent(in)
         490             :                               rho_ds_zm, invrs_rho_ds_zt,             & ! intent(in)
         491             :                               l_implemented, l_imp_sfc_momentum_flux, & ! intent(in)
         492           0 :                               lhs )                                     ! intent(out)
         493             :
         494             :       ! Decompose and back substitute for um and vm
         495           0 :       nrhs = 2
         496             :       call windm_edsclrm_solve( nz, ngrdcol, gr, nrhs, stats_metadata%iwindm_matrix_condt_num, & ! intent(in)
         497             :                                 tridiag_solve_method,                           & ! intent(in)
         498             :                                 stats_metadata,                                 & ! intent(in)
         499             :                                 stats_sfc, &                                      ! intent(inout)
         500             :                                 lhs, rhs, &                                       ! intent(inout)
         501           0 :                                 solution )                                        ! intent(out)
         502             :
         503             :       ! Check for singular matrices and bad LAPACK arguments
         504           0 :       if ( clubb_at_least_debug_level( 0 ) ) then
         505           0 :         if ( err_code == clubb_fatal_error ) then
         506           0 :           write(fstderr,*) "Fatal error solving for um/vm"
         507           0 :           return
         508             :         end if
         509             :       end if
         510             :
         511             :       !----------------------------------------------------------------
         512             :       ! Update zonal (west-to-east) component of mean wind, um
         513             :       !----------------------------------------------------------------
         514             :       !$acc parallel loop gang vector collapse(2) default(present)
         515           0 :       do k = 1, nz
         516           0 :         do i = 1, ngrdcol
         517           0 :           um(i,k) = solution(i,k,windm_edsclrm_um)
         518             :         end do
         519             :       end do
         520             :       !$acc end parallel loop
         521             :
         522             :       !----------------------------------------------------------------
         523             :       ! Update meridional (south-to-north) component of mean wind, vm
         524             :       !----------------------------------------------------------------
         525             :       !$acc parallel loop gang vector collapse(2) default(present)
         526           0 :       do k = 1, nz
         527           0 :         do i = 1, ngrdcol
         528           0 :           vm(i,k) = solution(i,k,windm_edsclrm_vm)
         529             :         end do
         530             :       end do
         531             :       !$acc end parallel loop
         532             :
         533           0 :       if ( stats_metadata%l_stats_samp ) then
         534             :
         535             :         !$acc update host( um, lhs_diff, lhs_ma_zt, invrs_rho_ds_zt, u_star_sqd, &
         536             :         !$acc              rho_ds_zm, wind_speed, vm )
         537             :
         538           0 :         do i = 1, ngrdcol
         539             :           ! Implicit contributions to um and vm
         540           0 :           call windm_edsclrm_implicit_stats( nz, windm_edsclrm_um, um(i,:),       & ! intent(in)
         541           0 :                                              gr%invrs_dzt(i,:),                   & ! intent(in)
         542           0 :                                              lhs_diff(:,i,:), lhs_ma_zt(:,i,:),   & ! intent(in)
         543           0 :                                              invrs_rho_ds_zt(i,:), u_star_sqd(i), & ! intent(in)
         544             :                                              rho_ds_zm(i,:), wind_speed(i,:),     & ! intent(in)
         545             :                                              l_imp_sfc_momentum_flux,             & ! intent(in)
         546             :                                              stats_metadata,                      & ! intent(in)
         547           0 :                                              stats_zt(i) )                          ! intent(inout)
         548             :
         549             :           call windm_edsclrm_implicit_stats( nz, windm_edsclrm_vm, vm(i,:),       & ! intent(in)
         550           0 :                                              gr%invrs_dzt(i,:),                   & ! intent(in)
         551             :                                              lhs_diff(:,i,:), lhs_ma_zt(:,i,:),   & ! intent(in)
         552             :                                              invrs_rho_ds_zt(i,:), u_star_sqd(i), & ! intent(in)
         553             :                                              rho_ds_zm(i,:), wind_speed(i,:),     & ! intent(in)
         554             :                                              l_imp_sfc_momentum_flux,             & ! intent(in)
         555             :                                              stats_metadata,                      & ! intent(in)
         556           0 :                                              stats_zt(i) )                          ! intent(inout)
         557             :         end do
         558             :       end if ! stats_metadata%l_stats_samp
         559             :
         560           0 :       if ( l_lmm_stepping ) then
         561             :         !$acc parallel loop gang vector collapse(2) default(present)
         562           0 :         do k = 1, nz
         563           0 :           do i = 1, ngrdcol
         564           0 :             um(i,k) = one_half * ( um_old(i,k) + um(i,k) )
         565           0 :             vm(i,k) = one_half * ( vm_old(i,k) + vm(i,k) )
         566             :           end do
         567             :         end do
         568             :         !$acc end parallel loop
         569             :       endif ! l_lmm_stepping ) then
         570             :
         571           0 :       if ( uv_sponge_damp_settings%l_sponge_damping ) then
         572             :
         573             :         !$acc update host( um, vm, um_ref, vm_ref )
         574             :
         575             :         ! _sponge_damp_settings and _sponge_damp_profile could potentially vary
         576             :         ! from column to column, but there is no column index in those variables.
         577             :         ! Thus this code is potentially unsafe when implemented in a host model,
         578             :         ! which is indicated by l_implemented = T
         579           0 :         if ( l_implemented ) then
         580             :           write(fstderr,*) "l_sponge_damping = T and l_implemented = T &
         581           0 :                             -- this is likely unsafe and considered fatal"
         582           0 :           err_code = clubb_fatal_error
         583           0 :           return
         584             :         end if
         585             :
         586           0 :         if ( stats_metadata%l_stats_samp ) then
         587           0 :           do i = 1, ngrdcol
         588           0 :             tmp_in(1) = 0.0_core_rknd
         589           0 :             tmp_in(2:nz) = um(i,2:nz)
         590             :             call stat_begin_update( nz, stats_metadata%ium_sdmp, tmp_in / dt, & ! intent(in)
         591           0 :                                     stats_zt(i) )           ! intent(inout)
         592           0 :             tmp_in(2:nz) = vm(i,2:nz)
         593             :             call stat_begin_update( nz, stats_metadata%ivm_sdmp, tmp_in / dt, & ! intent(in)
         594           0 :                                     stats_zt(i) )           ! intent(inout)
         595             :           end do
         596             :         end if
         597             :
         598           0 :         do i = 1, ngrdcol
         599           0 :           um(i,1:nz) = sponge_damp_xm( nz, dt, gr%zt(i,:), gr%zm(i,:), &
         600           0 :                                        um_ref(i,1:nz), um(i,1:nz), uv_sponge_damp_profile )
         601             :         end do
         602             :
         603           0 :         do i = 1, ngrdcol
         604           0 :           vm(i,1:nz) = sponge_damp_xm( nz, dt, gr%zt(i,:), gr%zm(i,:), &
         605           0 :                                        vm_ref(i,1:nz), vm(i,1:nz), uv_sponge_damp_profile )
         606             :         end do
         607             :
         608           0 :         if ( stats_metadata%l_stats_samp ) then
         609           0 :           do i = 1, ngrdcol
         610           0 :             tmp_in(1) = 0.0_core_rknd
         611           0 :             tmp_in(2:nz) = um(i,2:nz)
         612             :             call stat_end_update( nz, stats_metadata%ium_sdmp, tmp_in / dt, & ! intent(in)
         613           0 :                                   stats_zt(i) )           ! intent(inout)
         614           0 :             tmp_in(2:nz) = vm(i,2:nz)
         615             :             call stat_end_update( nz, stats_metadata%ivm_sdmp, tmp_in / dt, & ! intent(in)
         616           0 :                                   stats_zt(i) )           ! intent(inout)
         617             :           end do
         618             :         end if
         619             :
         620             :         !$acc update device( um, vm )
         621             :
         622             :       end if ! uv_sponge_damp_settings%l_sponge_damping
         623             :
         624             :       ! Second part of momentum (implicit component)
         625             :
         626             :       ! Solve for x'w' at all intermediate model levels.
         627             :       ! A Crank-Nicholson timestep is used.
         628             :       call calc_xpwp( nz, ngrdcol, gr, &
         629             :                       Km_zm_p_nu10, um, &
         630           0 :                       xpwp )
         631             :
         632             :       !$acc parallel loop gang vector collapse(2) default(present)
         633           0 :       do k = 2, nz-1
         634           0 :         do i = 1, ngrdcol
         635           0 :           upwp(i,k) = upwp(i,k) - one_half * xpwp(i,k)
         636             :         end do
         637             :       end do
         638             :       !$acc end parallel loop
         639             :
         640             :       call calc_xpwp( nz, ngrdcol, gr, &
         641             :                       Km_zm_p_nu10, vm, &
         642           0 :                       xpwp )
         643             :
         644             :       !$acc parallel loop gang vector collapse(2) default(present)
         645           0 :       do k = 2, nz-1
         646           0 :         do i = 1, ngrdcol
         647           0 :           vpwp(i,k) = vpwp(i,k) - one_half * xpwp(i,k)
         648             :         end do
         649             :       end do
         650             :       !$acc end parallel loop
         651             :
         652             :       ! Adjust um and vm if nudging is turned on.
         653           0 :       if ( l_uv_nudge ) then
         654             :
         655             :         ! Reflect nudging in budget
         656           0 :         if ( stats_metadata%l_stats_samp ) then
         657             :           !$acc update host( um, vm )
         658           0 :           do i = 1, ngrdcol
         659           0 :             tmp_in(1) = 0.0_core_rknd
         660           0 :             tmp_in(2:nz) = um(i,2:nz)
         661             :             call stat_begin_update( nz, stats_metadata%ium_ndg, tmp_in / dt, & ! intent(in)
         662           0 :                                     stats_zt(i) )          ! intent(inout)
         663           0 :             tmp_in(2:nz) = vm(i,2:nz)
         664             :             call stat_begin_update( nz, stats_metadata%ivm_ndg, tmp_in / dt, & ! intent(in)
         665           0 :                                     stats_zt(i) )          ! intent(inout)
         666             :           end do
         667             :         end if
         668             :
         669             :         !$acc parallel loop gang vector collapse(2) default(present)
         670           0 :         do k = 1, nz
         671           0 :           do i = 1, ngrdcol
         672           0 :             um(i,k) = um(i,k) - ( ( um(i,k) - um_ref(i,k) ) * (dt/ts_nudge) )
         673           0 :             vm(i,k) = vm(i,k) - ( ( vm(i,k) - vm_ref(i,k) ) * (dt/ts_nudge) )
         674             :           end do
         675             :         end do
         676             :         !$acc end parallel loop
         677             :
         678           0 :         if ( stats_metadata%l_stats_samp ) then
         679             :           !$acc update host( um, vm )
         680           0 :           do i = 1, ngrdcol
         681           0 :             tmp_in(1) = 0.0_core_rknd
         682           0 :             tmp_in(2:nz) = um(i,2:nz)
         683             :             call stat_end_update( nz, stats_metadata%ium_ndg, tmp_in / dt, & ! intent(in)
         684           0 :                                   stats_zt(i) )          ! intent(inout)
         685           0 :             tmp_in(2:nz) = vm(i,2:nz)
         686             :             call stat_end_update( nz, stats_metadata%ivm_ndg, tmp_in / dt, & ! intent(in)
         687           0 :                                   stats_zt(i) )          ! intent(inout)
         688             :           end do
         689             :         end if
         690             :
         691             :       end if ! l_uv_nudge
         692             :
         693           0 :       if ( stats_metadata%l_stats_samp ) then
         694             :         !$acc update host( um_ref, vm_ref )
         695           0 :         do i = 1, ngrdcol
         696           0 :           tmp_in(1) = 0.0_core_rknd
         697           0 :           tmp_in(2:nz) = um_ref(i,2:nz)
         698             :           call stat_update_var( stats_metadata%ium_ref, tmp_in, & ! intent(in)
         699           0 :                                 stats_zt(i) )         ! intent(inout)
         700           0 :           tmp_in(2:nz) = vm_ref(i,2:nz)
         701             :           call stat_update_var( stats_metadata%ivm_ref, tmp_in, & ! intent(in)
         702           0 :                                 stats_zt(i) )         ! intent(inout)
         703             :         end do
         704             :       end if
         705             :
         706             :       if ( order_windm < order_wp2_wp3 &
         707           0 :           .and. order_windm < order_xp2_xpyp ) then
         708           0 :         l_first_clip_ts = .true.
         709           0 :         l_last_clip_ts = .false.
         710             :       elseif ( order_windm > order_wp2_wp3 &
         711           0 :               .and. order_windm > order_xp2_xpyp ) then
         712           0 :         l_first_clip_ts = .false.
         713           0 :         l_last_clip_ts = .true.
         714             :       else
         715           0 :         l_first_clip_ts = .false.
         716           0 :         l_last_clip_ts = .false.
         717             :       endif
         718             :
         719           0 :       if ( l_tke_aniso ) then
         720             :
         721             :         ! Clipping for u'w'
         722             :         !
         723             :         ! Clipping u'w' at each vertical level, based on the
         724             :         ! correlation of u and w at each vertical level, such that:
         725             :         ! corr_(u,w) = u'w' / [ sqrt(u'^2) * sqrt(w'^2) ];
         726             :         ! -1 <= corr_(u,w) <= 1.
         727             :         !
         728             :         ! Since u'^2, w'^2, and u'w' are each advanced in different
         729             :         ! subroutines from each other in advance_clubb_core, clipping for u'w'
         730             :         ! has to be done three times during each timestep (once after each
         731             :         ! variable has been updated).
         732             :         ! This is the third instance of u'w' clipping.
         733             :         !l_first_clip_ts = .false.
         734             :         !l_last_clip_ts = .true.
         735             :         call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
         736             :                          l_last_clip_ts, dt, wp2, up2,                & ! intent(in)
         737             :                          l_predict_upwp_vpwp,                         & ! intent(in)
         738             :                          stats_metadata,                              & ! intent(in)
         739             :                          stats_zm,                                    & ! intent(inout)
         740           0 :                          upwp, upwp_chnge )                             ! intent(inout)
         741             :
         742             :         ! Clipping for v'w'
         743             :         !
         744             :         ! Clipping v'w' at each vertical level, based on the
         745             :         ! correlation of v and w at each vertical level, such that:
         746             :         ! corr_(v,w) = v'w' / [ sqrt(v'^2) * sqrt(w'^2) ];
         747             :         ! -1 <= corr_(v,w) <= 1.
         748             :         !
         749             :         ! Since v'^2, w'^2, and v'w' are each advanced in different
         750             :         ! subroutines from each other in advance_clubb_core, clipping for v'w'
         751             :         ! has to be done three times during each timestep (once after each
         752             :         ! variable has been updated).
         753             :         ! This is the third instance of v'w' clipping.
         754             :         !l_first_clip_ts = .false.
         755             :         !l_last_clip_ts = .true.
         756             :         call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
         757             :                          l_last_clip_ts, dt, wp2, vp2,                & ! intent(in)
         758             :                          l_predict_upwp_vpwp,                         & ! intent(in)
         759             :                          stats_metadata,                              & ! intent(in)
         760             :                          stats_zm,                                    & ! intent(inout)
         761           0 :                          vpwp, vpwp_chnge )                             ! intent(inout)
         762             :       else
         763             :
         764             :         ! intent(in) this case, it is assumed that
         765             :         !   u'^2 == v'^2 == w'^2, and the variables `up2' and `vp2' do not
         766             :         ! interact with any other variables.
         767             :         !l_first_clip_ts = .false.
         768             :         !l_last_clip_ts = .true.
         769             :         call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
         770             :                          l_last_clip_ts, dt, wp2, wp2,                & ! intent(in)
         771             :                          l_predict_upwp_vpwp,                         & ! intent(in)
         772             :                          stats_metadata,                              & ! intent(in)
         773             :                          stats_zm,                                    & ! intent(inout)
         774           0 :                          upwp, upwp_chnge )                             ! intent(inout)
         775             :
         776             :         call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
         777             :                          l_last_clip_ts, dt, wp2, wp2,                & ! intent(in)
         778             :                          l_predict_upwp_vpwp,                         & ! intent(in)
         779             :                          stats_metadata,                              & ! intent(in)
         780             :                          stats_zm,                                    & ! intent(inout)
         781           0 :                          vpwp, vpwp_chnge )                             ! intent(inout)
         782             :       endif ! l_tke_aniso
         783             :
         784             :       ! The values of um(1) and vm(1) are located below the model surface and
         785             :       ! do not affect the rest of the model.  The values of um(1) or vm(1) are
         786             :       ! simply set to the values of um(2) and vm(2), respectively, after the
         787             :       ! equation matrices has been solved.  Even though um and vm would sharply
         788             :       ! decrease to a value of 0 at the surface, this is done to avoid
         789             :       ! confusion on plots of the vertical profiles of um and vm.
         790             :       !$acc parallel loop gang vector default(present)
         791           0 :       do i = 1, ngrdcol
         792           0 :         um(i,1) = um(i,2)
         793           0 :         vm(i,1) = vm(i,2)
         794             :       end do
         795             :       !$acc end parallel loop
         796             :
         797             :     endif ! .not. l_predict_upwp_vpwp
         798             :
         799             :
         800     8935056 :     if ( l_perturbed_wind ) then
         801             :
         802             :       !----------------------------------------------------------------
         803             :       ! Prepare tridiagonal system for horizontal winds, um and vm
         804             :       !----------------------------------------------------------------
         805             :
         806             :       ! Momentum surface fluxes, u'w'|_sfc and v'w'|_sfc, are applied through
         807             :       ! an implicit method, such that:
         808             :       !    x'w'|_sfc = - ( u_star(t)^2 / wind_speed(t) ) * xm(t+1).
         809           0 :       l_imp_sfc_momentum_flux = .true.
         810             :
         811             :       ! Compute wind speed (use threshold "eps" to prevent divide-by-zero
         812             :       ! error).
         813             :       !$acc parallel loop gang vector collapse(2) default(present)
         814           0 :       do k = 1, nz
         815           0 :         do i = 1, ngrdcol
         816           0 :           wind_speed_pert(i,k) = max( sqrt( (um_pert(i,k))**2 + (vm_pert(i,k))**2 ), eps )
         817             :         end do
         818             :       end do
         819             :       !$acc end parallel loop
         820             :
         821             :       ! Compute u_star_sqd according to the definition of u_star.
         822             :       !$acc parallel loop gang vector default(present)
         823           0 :       do i = 1, ngrdcol
         824           0 :         u_star_sqd_pert(i) = sqrt( upwp_pert(i,1)**2 + vpwp_pert(i,1)**2 )
         825             :       end do
         826             :       !$acc end parallel loop
         827             :
         828             :       ! Compute the explicit portion of the um equation.
         829             :       ! Build the right-hand side vector.
         830             :       call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_um, dt,    & ! intent(in)
         831             :                               lhs_diff, um_pert, um_tndcy,              & ! intent(in)
         832             :                               rho_ds_zm, invrs_rho_ds_zt,               & ! intent(in)
         833             :                               l_imp_sfc_momentum_flux, upwp_pert(:,1),  & ! intent(in)
         834             :                               stats_metadata,                           & ! intent(in)
         835             :                               stats_zt,                                 & ! intent(inout)
         836           0 :                               rhs(:,:,windm_edsclrm_um) )                 ! intent(out)
         837             :
         838             :       ! Compute the explicit portion of the vm equation.
         839             :       ! Build the right-hand side vector.
         840             :       call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_vm, dt,    & ! intent(in)
         841             :                               lhs_diff, vm_pert, vm_tndcy,              & ! intent(in)
         842             :                               rho_ds_zm, invrs_rho_ds_zt,               & ! intent(in)
         843             :                               l_imp_sfc_momentum_flux, vpwp_pert(:,1),  & ! intent(in)
         844             :                               stats_metadata,                           & ! intent(in)
         845             :                               stats_zt,                                 & ! intent(inout)
         846           0 :                               rhs(:,:,windm_edsclrm_vm) )                 ! intent(out)
         847             :
         848             :       ! Store momentum flux (explicit component)
         849             :
         850             :       ! The surface flux, x'w'(1) = x'w'|_sfc, is set elsewhere in the model.
         851             :       !      upwp(1) = upwp_sfc
         852             :       !      vpwp(1) = vpwp_sfc
         853             :
         854             :       ! Solve for x'w' at all intermediate model levels.
         855             :       ! A Crank-Nicholson timestep is used.
         856             :       call calc_xpwp( nz, ngrdcol, gr, &
         857             :                       Km_zm_p_nu10, um_pert, &
         858           0 :                       xpwp )
         859             :
         860             :       !$acc parallel loop gang vector collapse(2) default(present)
         861           0 :       do k = 2, nz-1
         862           0 :         do i = 1, ngrdcol
         863           0 :           upwp_pert(i,k) = -one_half * xpwp(i,k)
         864             :         end do
         865             :       end do
         866             :       !$acc end parallel loop
         867             :
         868             :       call calc_xpwp( nz, ngrdcol, gr, &
         869             :                       Km_zm_p_nu10, vm_pert, &
         870           0 :                       xpwp )
         871             :
         872             :       !$acc parallel loop gang vector collapse(2) default(present)
         873           0 :       do k = 2, nz-1
         874           0 :         do i = 1, ngrdcol
         875           0 :           vpwp_pert(i,k) = -one_half * xpwp(i,k)
         876             :         end do
         877             :       end do
         878             :       !$acc end parallel loop
         879             :
         880             :       ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
         881             :       ! means that x'w' at the top model level is 0,
         882             :       ! since x'w' = - K_zm * d(xm)/dz.
         883             :       !$acc parallel loop gang vector default(present)
         884           0 :       do i = 1, ngrdcol
         885           0 :         upwp_pert(i,nz) = zero
         886           0 :         vpwp_pert(i,nz) = zero
         887             :       end do
         888             :       !$acc end parallel loop
         889             :
         890             :       ! Compute the implicit portion of the um and vm equations.
         891             :       ! Build the left-hand side matrix.
         892             :       call windm_edsclrm_lhs( nz, ngrdcol, gr, dt,                        & ! intent(in)
         893             :                               lhs_ma_zt, lhs_diff,                        & ! intent(in)
         894             :                               wind_speed_pert, u_star_sqd_pert,           & ! intent(in)
         895             :                               rho_ds_zm, invrs_rho_ds_zt,                 & ! intent(in)
         896             :                               l_implemented, l_imp_sfc_momentum_flux,     & ! intent(in)
         897           0 :                               lhs )                                         ! intent(out)
         898             :
         899             :       ! Decompose and back substitute for um and vm
         900           0 :       nrhs = 2
         901             :       call windm_edsclrm_solve( nz, ngrdcol, gr, nrhs, stats_metadata%iwindm_matrix_condt_num, & ! intent(in)
         902             :                                 tridiag_solve_method,                           & ! intent(in)
         903             :                                 stats_metadata,                                 & ! intent(in)
         904             :                                 stats_sfc, &                                      ! intent(in)
         905             :                                 lhs, rhs, &                                       ! intent(inout)
         906           0 :                                 solution )                                        ! intent(out)
         907             :
         908             :       ! Check for singular matrices and bad LAPACK arguments
         909           0 :       if ( clubb_at_least_debug_level( 0 ) ) then
         910           0 :         if ( err_code == clubb_fatal_error ) then
         911           0 :           write(fstderr,*) "Fatal error solving for um_pert/vm_pert"
         912           0 :           return
         913             :         endif
         914             :       endif
         915             :
         916             :       !----------------------------------------------------------------
         917             :       ! Update zonal (west-to-east) component of mean wind, um
         918             :       !----------------------------------------------------------------
         919             :       !$acc parallel loop gang vector collapse(2) default(present)
         920           0 :       do k = 1, nz
         921           0 :         do i = 1, ngrdcol
         922           0 :           um_pert(i,k) = solution(i,k,windm_edsclrm_um)
         923             :         end do
         924             :       end do
         925             :       !$acc end parallel loop
         926             :
         927             :       !----------------------------------------------------------------
         928             :       ! Update meridional (south-to-north) component of mean wind, vm
         929             :       !----------------------------------------------------------------
         930             :       !$acc parallel loop gang vector collapse(2) default(present)
         931           0 :       do k = 1, nz
         932           0 :         do i = 1, ngrdcol
         933           0 :           vm_pert(i,k) = solution(i,k,windm_edsclrm_vm)
         934             :         end do
         935             :       end do
         936             :       !$acc end parallel loop
         937             :
         938             :       ! Second part of momentum (implicit component)
         939             :
         940             :       ! Solve for x'w' at all intermediate model levels.
         941             :       ! A Crank-Nicholson timestep is used.
         942             :       call calc_xpwp( nz, ngrdcol, gr, &
         943             :                       Km_zm_p_nu10, um_pert, &
         944           0 :                       xpwp )
         945             :
         946             :       !$acc parallel loop gang vector collapse(2) default(present)
         947           0 :       do k = 2, nz-1
         948           0 :         do i = 1, ngrdcol
         949           0 :           upwp_pert(i,k) = upwp_pert(i,k) - one_half * xpwp(i,k)
         950             :         end do
         951             :       end do
         952             :       !$acc end parallel loop
         953             :
         954             :       call calc_xpwp( nz, ngrdcol, gr, &
         955             :                       Km_zm_p_nu10, vm_pert, &
         956           0 :                       xpwp )
         957             :
         958             :       !$acc parallel loop gang vector collapse(2) default(present)
         959           0 :       do k = 2, nz-1
         960           0 :         do i = 1, ngrdcol
         961           0 :           vpwp_pert(i,k) = vpwp_pert(i,k) - one_half * xpwp(i,k)
         962             :         end do
         963             :       end do
         964             :       !$acc end parallel loop
         965             :
         966           0 :       if ( l_tke_aniso ) then
         967             :
         968             :         ! Clipping for u'w'
         969             :         !
         970             :         ! Clipping u'w' at each vertical level, based on the
         971             :         ! correlation of u and w at each vertical level, such that:
         972             :         ! corr_(u,w) = u'w' / [ sqrt(u'^2) * sqrt(w'^2) ];
         973             :         ! -1 <= corr_(u,w) <= 1.
         974             :         !
         975             :         ! Since u'^2, w'^2, and u'w' are each advanced in different
         976             :         ! subroutines from each other in advance_clubb_core, clipping for u'w'
         977             :         ! has to be done three times during each timestep (once after each
         978             :         ! variable has been updated).
         979             :         ! This is the third instance of u'w' clipping.
         980           0 :         l_first_clip_ts = .false.
         981           0 :         l_last_clip_ts = .true.
         982             :         call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
         983             :                          l_last_clip_ts, dt, wp2, up2,                & ! intent(in)
         984             :                          l_predict_upwp_vpwp,                         & ! intent(in)
         985             :                          stats_metadata,                              & ! intent(in)
         986             :                          stats_zm,                                    & ! intent(inout)
         987           0 :                          upwp_pert, upwp_chnge )                        ! intent(inout)
         988             :
         989             :         ! Clipping for v'w'
         990             :         !
         991             :         ! Clipping v'w' at each vertical level, based on the
         992             :         ! correlation of v and w at each vertical level, such that:
         993             :         ! corr_(v,w) = v'w' / [ sqrt(v'^2) * sqrt(w'^2) ];
         994             :         ! -1 <= corr_(v,w) <= 1.
         995             :         !
         996             :         ! Since v'^2, w'^2, and v'w' are each advanced in different
         997             :         ! subroutines from each other in advance_clubb_core, clipping for v'w'
         998             :         ! has to be done three times during each timestep (once after each
         999             :         ! variable has been updated).
        1000             :         ! This is the third instance of v'w' clipping.
        1001           0 :         l_first_clip_ts = .false.
        1002           0 :         l_last_clip_ts = .true.
        1003             :         call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
        1004             :                          l_last_clip_ts, dt, wp2, vp2,                & ! intent(in)
        1005             :                          l_predict_upwp_vpwp,                         & ! intent(in)
        1006             :                          stats_metadata,                              & ! intent(in)
        1007             :                          stats_zm,                                    & ! intent(inout)
        1008           0 :                          vpwp_pert, vpwp_chnge )                        ! intent(inout)
        1009             :       else
        1010             :
        1011             :         ! intent(in) this case, it is assumed that
        1012             :         !   u'^2 == v'^2 == w'^2, and the variables `up2' and `vp2' do not
        1013             :         ! interact with any other variables.
        1014           0 :         l_first_clip_ts = .false.
        1015           0 :         l_last_clip_ts = .true.
        1016             :         call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
        1017             :                          l_last_clip_ts, dt, wp2, wp2,                & ! intent(in)
        1018             :                          l_predict_upwp_vpwp,                         & ! intent(in)
        1019             :                          stats_metadata,                              & ! intent(in)
        1020             :                          stats_zm,                                    & ! intent(inout)
        1021           0 :                          upwp_pert, upwp_chnge )                        ! intent(inout)
        1022             :
        1023             :         call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
        1024             :                          l_last_clip_ts, dt, wp2, wp2,                & ! intent(in)
        1025             :                          l_predict_upwp_vpwp,                         & ! intent(in)
        1026             :                          stats_metadata,                              & ! intent(in)
        1027             :                          stats_zm,                                    & ! intent(inout)
        1028           0 :                          vpwp_pert, vpwp_chnge )                        ! intent(inout)
        1029             :
        1030             :       end if ! l_tke_aniso
        1031             :
        1032             :       ! The values of um(1) and vm(1) are located below the model surface and
        1033             :       ! do not affect the rest of the model.  The values of um(1) or vm(1) are
        1034             :       ! simply set to the values of um(2) and vm(2), respectively, after the
        1035             :       ! equation matrices has been solved.  Even though um and vm would sharply
        1036             :       ! decrease to a value of 0 at the surface, this is done to avoid
        1037             :       ! confusion on plots of the vertical profiles of um and vm.
        1038             :       !$acc parallel loop gang vector default(present)
        1039           0 :       do i = 1, ngrdcol
        1040           0 :         um_pert(i,1) = um_pert(i,2)
        1041           0 :         vm_pert(i,1) = vm_pert(i,2)
        1042             :       end do
        1043             :       !$acc end parallel loop
        1044             :
        1045             :     end if ! l_perturbed_wind
        1046             :
        1047             :     !----------------------------------------------------------------
        1048             :     ! Prepare tridiagonal system for eddy-scalars
        1049             :     !----------------------------------------------------------------
        1050             :
        1051     8935056 :     if ( edsclr_dim > 0 ) then
        1052             :
        1053     8935056 :       Kmh_zt(:,:) = zm2zt( nz, ngrdcol, gr, Kmh_zm(:,:), zero )
        1054             :
        1055             :       ! Calculate diffusion terms
        1056             :       call diffusion_zt_lhs( nz, ngrdcol, gr, Kmh_zm, Kmh_zt, nu_zero,  & ! intent(in)
        1057             :                              invrs_rho_ds_zt, rho_ds_zm,                & ! intent(in)
        1058     8935056 :                              lhs_diff )                                   ! intent(out)
        1059             :
        1060             :       ! The lower boundary condition needs to be applied here at level 2.
        1061             :       if ( .not. l_upwind_Kh_dp_term ) then
        1062             :
        1063             :         !$acc parallel loop gang vector default(present)
        1064   149194656 :         do i = 1, ngrdcol
        1065             :           ! Thermodynamic superdiagonal: [ x xm(k+1,<t+1>) ]
        1066   280519200 :           lhs_diff(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
        1067             :                                     * ( Kmh_zm(i,2) + nu_zero(i) ) &
        1068   280519200 :                                     * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
        1069             :
        1070             :           ! Thermodynamic main diagonal: [ x xm(k,<t+1>) ]
        1071           0 :           lhs_diff(k_tdiag,i,2) = + gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
        1072             :                                   * ( Kmh_zm(i,2) + nu_zero(i) ) &
        1073   140259600 :                                   * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
        1074             :
        1075             :           ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
        1076   149194656 :           lhs_diff(km1_tdiag,i,2) = zero
        1077             :         end do
        1078             :         !$acc end parallel loop
        1079             :
        1080             :       else
        1081             :
        1082             :         !$acc parallel loop gang vector default(present)
        1083             :         do i = 1, ngrdcol
        1084             :           ! Thermodynamic superdiagonal: [ x xm(k+1,<t+1>) ]
        1085             :           lhs_diff(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) &
        1086             :                                       * ( Kmh_zt(i,2) + nu_zero(i) ) * gr%invrs_dzm(i,2)
        1087             :
        1088             :           ! Thermodynamic main diagonal: [ x xm(k,<t+1>) ]
        1089             :           lhs_diff(k_tdiag,i,2) = + gr%invrs_dzt(i,2) &
        1090             :                                     * ( Kmh_zt(i,2) + nu_zero(i) ) * gr%invrs_dzm(i,2)
        1091             :
        1092             :           ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
        1093             :           lhs_diff(km1_tdiag,i,2) = zero
        1094             :         end do
        1095             :         !$acc end parallel loop
        1096             :
        1097             :       end if
        1098             :
        1099     8935056 :       if ( l_lmm_stepping ) then
        1100             :         !$acc parallel loop gang vector collapse(3) default(present)
        1101           0 :         do j = 1, edsclr_dim
        1102           0 :           do k = 1, nz
        1103           0 :             do i = 1, ngrdcol
        1104           0 :               edsclrm_old(i,k,j) = edsclrm(i,k,j)
        1105             :             end do
        1106             :           end do
        1107             :         end do
        1108             :         !$acc end parallel loop
        1109             :       endif ! l_lmm_stepping
        1110             :
        1111             :       ! Eddy-scalar surface fluxes, x'w'|_sfc, are applied through an explicit
        1112             :       ! method.
        1113     8935056 :       l_imp_sfc_momentum_flux = .false.
        1114             :
        1115             :       ! Compute the explicit portion of eddy scalar equation.
        1116             :       ! Build the right-hand side vector.
        1117             :       ! Because of statistics, we have to use a DO rather than a FORALL here
        1118             :       ! -dschanen 7 Oct 2008
        1119   125090784 :       do edsclr = 1, edsclr_dim
        1120             :         call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_scalar, dt,      & ! intent(in)
        1121             :                                 lhs_diff, edsclrm(:,:,edsclr),                  & ! intent(in)
        1122             :                                 edsclrm_forcing(:,:,edsclr),                    & ! intent(in)
        1123             :                                 rho_ds_zm, invrs_rho_ds_zt,                     & ! intent(in)
        1124             :                                 l_imp_sfc_momentum_flux, wpedsclrp(:,1,edsclr), & ! intent(in)
        1125             :                                 stats_metadata,                                 & ! intent(in)
        1126             :                                 stats_zt,                                       & ! intent(inout)
        1127   125090784 :                                 rhs(:,:,edsclr) )                                 ! intent(out)
        1128             :       enddo
        1129             :
        1130             :
        1131             :       ! Store momentum flux (explicit component)
        1132             :
        1133             :       ! The surface flux, x'w'(1) = x'w'|_sfc, is set elsewhere in the model.
        1134             :       ! wpedsclrp(1,1:edsclr_dim) =  wpedsclrp_sfc(1:edsclr_dim)
        1135             :
        1136             :       ! Solve for x'w' at all intermediate model levels.
        1137             :       ! A Crank-Nicholson timestep is used
        1138   125090784 :       do edsclr = 1, edsclr_dim
        1139             :
        1140             :         call calc_xpwp( nz, ngrdcol, gr, &
        1141             :                         Km_zm_p_nu10, edsclrm(:,:,edsclr), &
        1142   116155728 :                         xpwp )
        1143             :
        1144             :         !$acc parallel loop gang vector collapse(2) default(present)
        1145  9766016208 :         do k = 2, nz-1
        1146 >16109*10^7 :           do i = 1, ngrdcol
        1147 >16098*10^7 :             wpedsclrp(i,k,edsclr) = -one_half * xpwp(i,k)
        1148             :           end do
        1149             :         end do
        1150             :         !$acc end parallel loop
        1151             :       end do
        1152             :
        1153             :       ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
        1154             :       ! means that x'w' at the top model level is 0,
        1155             :       ! since x'w' = - K_zm * d(xm)/dz.
        1156             :       !$acc parallel loop gang vector collapse(2) default(present)
        1157   125090784 :       do j = 1, edsclr_dim
        1158  1948465584 :         do i = 1, ngrdcol
        1159  1939530528 :           wpedsclrp(i,nz,j) = zero
        1160             :         end do
        1161             :       end do
        1162             :       !$acc end parallel loop
        1163             :
        1164             :       ! Compute the implicit portion of the xm (eddy-scalar) equations.
        1165             :       ! Build the left-hand side matrix.
        1166             :       call windm_edsclrm_lhs( nz, ngrdcol, gr, dt,                    & ! intent(in)
        1167             :                               lhs_ma_zt, lhs_diff,                    & ! intent(in)
        1168             :                               wind_speed, u_star_sqd,                 & ! intent(in)
        1169             :                               rho_ds_zm, invrs_rho_ds_zt,             & ! intent(in)
        1170             :                               l_implemented, l_imp_sfc_momentum_flux, & ! intent(in)
        1171     8935056 :                               lhs )                                     ! intent(out)
        1172             :
        1173             :       ! Decompose and back substitute for all eddy-scalar variables
        1174             :       call windm_edsclrm_solve( nz, ngrdcol, gr, edsclr_dim, 0, & ! intent(in)
        1175             :                                 tridiag_solve_method,           & ! intent(in)
        1176             :                                 stats_metadata,                 & ! intent(in)
        1177             :                                 stats_sfc,                      & ! intent(inout)
        1178             :                                 lhs, rhs,                       & ! intent(inout)
        1179     8935056 :                                 solution )                        ! intent(out)
        1180             :
        1181     8935056 :       if ( clubb_at_least_debug_level( 0 ) ) then
        1182     8935056 :         if ( err_code == clubb_fatal_error ) then
        1183           0 :           write(fstderr,*) "Fatal error solving for eddsclrm"
        1184             :         end if
        1185             :       end if
        1186             :
        1187             :       !----------------------------------------------------------------
        1188             :       ! Update Eddy-diff. Passive Scalars
        1189             :       !----------------------------------------------------------------
        1190             :       !$acc parallel loop gang vector collapse(3) default(present)
        1191   125090784 :       do j = 1, edsclr_dim
        1192  9998327664 :         do k = 1, nz
        1193 >16497*10^7 :           do i = 1, ngrdcol
        1194 >16486*10^7 :             edsclrm(i,k,j) = solution(i,k,j)
        1195             :           end do
        1196             :         end do
        1197             :       end do
        1198             :       !$acc end parallel loop
        1199             :
        1200     8935056 :       if ( l_lmm_stepping ) then
        1201             :         !$acc parallel loop gang vector collapse(3) default(present)
        1202           0 :         do j = 1, edsclr_dim
        1203           0 :           do k = 1, nz
        1204           0 :             do i = 1, ngrdcol
        1205           0 :               edsclrm(i,k,j) = one_half * ( edsclrm_old(i,k,j) + edsclrm(i,k,j) )
        1206             :             end do
        1207             :           end do
        1208             :         end do
        1209             :         !$acc end parallel loop
        1210             :       endif ! l_lmm_stepping
        1211             :
        1212             :       ! Second part of momentum (implicit component)
        1213             :
        1214             :       ! Solve for x'w' at all intermediate model levels.
        1215             :       ! A Crank-Nicholson timestep is used.
        1216   125090784 :       do edsclr = 1, edsclr_dim
        1217             :
        1218             :         call calc_xpwp( nz, ngrdcol, gr, &
        1219             :                         Kmh_zm, edsclrm(:,:,edsclr), &
        1220   116155728 :                         xpwp )
        1221             :
        1222             :         !$acc parallel loop gang vector collapse(2) default(present)
        1223  9766016208 :         do k = 2, nz-1
        1224 >16109*10^7 :           do i = 1, ngrdcol
        1225 >16098*10^7 :             wpedsclrp(i,k,edsclr) = -one_half * xpwp(i,k)
        1226             :           end do
        1227             :         end do
        1228             :         !$acc end parallel loop
        1229             :       end do
        1230             :
        1231             :       ! Note that the w'edsclr' terms are not clipped, since we don't compute
        1232             :       ! the variance of edsclr anywhere. -dschanen 7 Oct 2008
        1233             :
        1234             :       ! The value of edsclrm(1) is located below the model surface and does not
        1235             :       ! effect the rest of the model.  The value of edsclrm(1) is simply set to
        1236             :       ! the value of edsclrm(2) after the equation matrix has been solved.
        1237             :       !$acc parallel loop gang vector collapse(2) default(present)
        1238   125090784 :       do edsclr = 1, edsclr_dim
        1239  1948465584 :         do i = 1, ngrdcol
        1240  1939530528 :           edsclrm(i,1,edsclr) = edsclrm(i,2,edsclr)
        1241             :         end do
        1242             :       end do
        1243             :       !$acc end parallel loop
        1244             :
        1245             :     endif
        1246             :
        1247     8935056 :     if ( clubb_at_least_debug_level( 0 ) ) then
        1248     8935056 :         if ( err_code == clubb_fatal_error ) then
        1249             :
        1250             :           !$acc update host( wm_zt, Km_zm, ug, vg, um_ref, vm_ref, wp2, &
        1251             :           !$acc              up2, vp2, um_forcing, vm_forcing, edsclrm_forcing, &
        1252             :           !$acc              fcor, um_old, um, vm_old, vm, edsclrm_old, &
        1253             :           !$acc              edsclrm, upwp, vpwp, wpedsclrp )
        1254             :
        1255           0 :           write(fstderr,*) "Error in advance_windm_edsclrm"
        1256             :
        1257           0 :           write(fstderr,*) "intent(in)"
        1258             :
        1259           0 :           write(fstderr,*) "dt = ", dt
        1260           0 :           write(fstderr,*) "wm_zt = ", wm_zt
        1261           0 :           write(fstderr,*) "Km_zm = ", Km_zm
        1262           0 :           write(fstderr,*) "ug = ", ug
        1263           0 :           write(fstderr,*) "vg = ", vg
        1264           0 :           write(fstderr,*) "um_ref = ", um_ref
        1265           0 :           write(fstderr,*) "vm_ref = ", vm_ref
        1266           0 :           write(fstderr,*) "wp2 = ", wp2
        1267           0 :           write(fstderr,*) "up2 = ", up2
        1268           0 :           write(fstderr,*) "vp2 = ", vp2
        1269           0 :           write(fstderr,*) "um_forcing = ", um_forcing
        1270           0 :           write(fstderr,*) "vm_forcing = ", vm_forcing
        1271           0 :           do edsclr = 1, edsclr_dim
        1272           0 :             write(fstderr,*) "edsclrm_forcing # = ",edsclr, edsclrm_forcing
        1273             :           end do
        1274           0 :           write(fstderr,*) "fcor = ", fcor
        1275           0 :           write(fstderr,*) "l_implemented = ", l_implemented
        1276             :
        1277           0 :           write(fstderr,*) "intent(inout)"
        1278             :
        1279           0 :           if ( l_lmm_stepping ) &
        1280           0 :              write(fstderr,*) "um (pre-solve) = ", um_old
        1281           0 :           write(fstderr,*) "um = ", um
        1282           0 :           if ( l_lmm_stepping ) &
        1283           0 :              write(fstderr,*) "vm (pre-solve) = ", vm_old
        1284           0 :           write(fstderr,*) "vm = ", vm
        1285           0 :           do edsclr = 1, edsclr_dim
        1286           0 :             if ( l_lmm_stepping ) &
        1287           0 :                write(fstderr,*) "edsclrm (pre-solve) # ", edsclr, "=", edsclrm_old(:,:,edsclr)
        1288           0 :             write(fstderr,*) "edsclrm # ", edsclr, "=", edsclrm(:,:,edsclr)
        1289             :           end do
        1290           0 :           write(fstderr,*) "upwp = ", upwp
        1291           0 :           write(fstderr,*) "vpwp = ", vpwp
        1292           0 :           write(fstderr,*) "wpedsclrp = ", wpedsclrp
        1293             :
        1294           0 :           return
        1295             :         end if
        1296             :     end if
        1297             :
        1298             :     !$acc exit data delete( um_old, vm_old, um_tndcy, vm_tndcy, &
        1299             :     !$acc                    upwp_chnge, vpwp_chnge, lhs, rhs, solution, wind_speed, &
        1300             :     !$acc                    wind_speed_pert, u_star_sqd, u_star_sqd_pert, &
        1301             :     !$acc                    nu_zero, lhs_diff, lhs_ma_zt, Km_zt, Kmh_zt, &
        1302             :     !$acc                    Km_zm_p_nu10, xpwp )
        1303             :
        1304             :     !$acc exit data if( edsclr_dim > 0) delete( edsclrm_old )
        1305             :
        1306             :     return
        1307             :
        1308             :   end subroutine advance_windm_edsclrm
        1309             :
        1310             :   !=============================================================================
        1311     8935056 :   subroutine windm_edsclrm_solve( nz, ngrdcol, gr, nrhs, ixm_matrix_condt_num, &
        1312             :                                   tridiag_solve_method, &
        1313             :                                   stats_metadata, &
        1314     8935056 :                                   stats_sfc, &
        1315     8935056 :                                   lhs, rhs, solution )
        1316             :
        1317             :     ! Note:  In the "Description" section of this subroutine, the variable
        1318             :     !        "invrs_dzm" will be written as simply "dzm", and the variable
        1319             :     !        "invrs_dzt" will be written as simply "dzt".  This is being done as
        1320             :     !        as device to save space and to make some parts of the description
        1321             :     !        more readable.  This change does not pertain to the actual code.
        1322             :
        1323             :     ! Description:
        1324             :     ! Solves the horizontal wind or eddy-scalar time-tendency equation, and
        1325             :     ! diagnoses the turbulent flux.  A Crank-Nicholson time-stepping algorithm
        1326             :     ! is used in solving the turbulent advection term and in diagnosing the
        1327             :     ! turbulent flux.
        1328             :     !
        1329             :     ! The rate of change of an eddy-scalar quantity, xm, is:
        1330             :     !
        1331             :     ! d(xm)/dt = - w * d(xm)/dz - (1/rho_ds) * d( rho_ds * x'w' )/dz
        1332             :     !            + xm_forcings.
        1333             :     !
        1334             :     !
        1335             :     ! The Turbulent Advection Term
        1336             :     ! ----------------------------
        1337             :     !
        1338             :     ! The above equation contains a turbulent advection term:
        1339             :     !
        1340             :     ! - (1/rho_ds) * d( rho_ds * x'w' )/dz;
        1341             :     !
        1342             :     ! where the momentum flux, x'w', is closed using a down gradient approach:
        1343             :     !
        1344             :     ! x'w' = - K_zm * d(xm)/dz.
        1345             :     !
        1346             :     ! The turbulent advection term becomes:
        1347             :     !
        1348             :     ! + (1/rho_ds) * d [ rho_ds * K_zm * d(xm)/dz ] / dz;
        1349             :     !
        1350             :     ! which is the same as a standard eddy-diffusion term (if "rho_ds * K_zm" in
        1351             :     ! the term above is substituted for "K_zm" in a standard eddy-diffusion
        1352             :     ! term, and if the standard eddy-diffusion term is multiplied by
        1353             :     ! "1/rho_ds").  Thus, the turbulent advection term is treated and solved in
        1354             :     ! the same way that a standard eddy-diffusion term would be solved.  The
        1355             :     ! term is discretized as follows:
        1356             :     !
        1357             :     ! The values of xm are found on the thermodynamic levels, while the values
        1358             :     ! of K_zm are found on the momentum levels.  Additionally, the values of
        1359             :     ! rho_ds_zm are found on the momentum levels, and the values of
        1360             :     ! invrs_rho_ds_zt are found on the thermodynamic levels.  The
        1361             :     ! derivatives (d/dz) of xm are taken over the intermediate momentum levels.
        1362             :     ! At the intermediate momentum levels, d(xm)/dz is multiplied by K_zm and by
        1363             :     ! rho_ds_zm.  Then, the derivative of the whole mathematical expression is
        1364             :     ! taken over the central thermodynamic level, where it is multiplied by
        1365             :     ! invrs_rho_ds_zt, which yields the desired result.
        1366             :     !
        1367             :     ! ---xm(kp1)----------------------------------------------------- t(k+1)
        1368             :     !
        1369             :     ! ===========d(xm)/dz===K_zm(k)=====rho_ds_zm(k)================= m(k)
        1370             :     !
        1371             :     ! ---xm(k)---invrs_rho_ds_zt---d[rho_ds_zm*K_zm*d(xm)/dz]/dz----- t(k)
        1372             :     !
        1373             :     ! ===========d(xm)/dz===K_zm(km1)===rho_ds_zm(km1)=============== m(k-1)
        1374             :     !
        1375             :     ! ---xm(km1)----------------------------------------------------- t(k-1)
        1376             :     !
        1377             :     ! The vertical indices t(k+1), m(k), t(k), m(k-1), and t(k-1) correspond
        1378             :     ! with altitudes zt(k+1), zm(k), zt(k), zm(k-1), and zt(k-1), respectively.
        1379             :     ! The letter "t" is used for thermodynamic levels and the letter "m" is used
        1380             :     ! for momentum levels.
        1381             :     !
        1382             :     ! dzt(k)   = 1 / ( zm(k) - zm(k-1) )
        1383             :     ! dzm(k)   = 1 / ( zt(k+1) - zt(k) )
        1384             :     ! dzm(k-1) = 1 / ( zt(k) - zt(k-1) )
        1385             :     !
        1386             :     ! The vertically discretized form of the turbulent advection term (treated
        1387             :     ! as an eddy diffusion term) is written out as:
        1388             :     !
        1389             :     ! + invrs_rho_ds_zt(k)
        1390             :     !   * dzt(k)
        1391             :     !     * [   rho_ds_zm(k) * K_zm(k) * dzm(k) * ( xm(k+1) - xm(k) )
        1392             :     !         - rho_ds_zm(k-1) * K_zm(k-1) * dzm(k-1) * ( xm(k) - xm(k-1) ) ].
        1393             :     !
        1394             :     ! For this equation, a Crank-Nicholson (semi-implicit) diffusion scheme is
        1395             :     ! used to solve the (1/rho_ds) * d [ rho_ds * K_zm * d(xm)/dz ] / dz
        1396             :     ! eddy-diffusion term.  The discretized implicit portion of the term is
        1397             :     ! written out as:
        1398             :     !
        1399             :     ! + (1/2) * invrs_rho_ds_zt(k)
        1400             :     !   * dzt(k)
        1401             :     !     * [   rho_ds_zm(k) * K_zm(k)
        1402             :     !           * dzm(k) * ( xm(k+1,<t+1>) - xm(k,<t+1>) )
        1403             :     !         - rho_ds_zm(k-1) * K_zm(k-1)
        1404             :     !           * dzm(k-1) * ( xm(k,<t+1>) - xm(k-1,<t+1>) ) ].
        1405             :     !
        1406             :     ! Note:  When the implicit term is brought over to the left-hand side,
        1407             :     !        the sign is reversed and the leading "+" in front of the term
        1408             :     !        is changed to a "-".
        1409             :     !
        1410             :     ! The discretized explicit portion of the term is written out as:
        1411             :     !
        1412             :     ! + (1/2) * invrs_rho_ds_zt(k)
        1413             :     !   * dzt(k)
        1414             :     !     * [   rho_ds_zm(k) * K_zm(k)
        1415             :     !           * dzm(k) * ( xm(k+1,<t>) - xm(k,<t>) )
        1416             :     !         - rho_ds_zm(k-1) * K_zm(k-1)
        1417             :     !           * dzm(k-1) * ( xm(k,<t>) - xm(k-1,<t>) ) ].
        1418             :     !
        1419             :     ! Timestep index (t) stands for the index of the current timestep, while
        1420             :     ! timestep index (t+1) stands for the index of the next timestep, which is
        1421             :     ! being advanced to in solving the d(xm)/dt equation.
        1422             :     !
        1423             :     !
        1424             :     ! Boundary Conditions:
        1425             :     !
        1426             :     ! An eddy-scalar quantity is not allowed to flux out the upper boundary.
        1427             :     ! Thus, a zero-flux boundary condition is used for the upper boundary in the
        1428             :     ! eddy-diffusion equation.
        1429             :     !
        1430             :     ! The lower boundary condition is much more complicated.  It is neither a
        1431             :     ! zero-flux nor a fixed-point boundary condition.  Rather, it is a
        1432             :     ! fixed-flux boundary condition.  This term is a turbulent advection term,
        1433             :     ! but with the eddy-scalars, the only value of x'w' relevant in solving the
        1434             :     ! d(xm)/dt equation is the value of x'w' at the surface (the first momentum
        1435             :     ! level), which is written as x'w'|_sfc.
        1436             :     !
        1437             :     ! 1) x'w' surface flux; generalized explicit form
        1438             :     !
        1439             :     !    The x'w' surface flux is applied to the d(xm)/dt equation through the
        1440             :     !    turbulent advection term, which is:
        1441             :     !
        1442             :     !    - (1/rho_ds) * d( rho_ds * x'w' )/dz.
        1443             :     !
        1444             :     !    At most vertical levels, a substitution can be made for x'w', such
        1445             :     !    that:
        1446             :     !
        1447             :     !    x'w' = - K_zm * d(xm)/dz.
        1448             :     !
        1449             :     !    However, the same substitution cannot be made at the surface (momentum
        1450             :     !    level 1), as x'w'|_sfc is a surface flux that is explicitly computed
        1451             :     !    elsewhere in the model code.
        1452             :     !
        1453             :     !    The lower boundary condition, which in this case needs to be applied to
        1454             :     !    the d(xm)/dt equation at level 2, is discretized as follows:
        1455             :     !
        1456             :     !    --xm(3)------------------------------------------------------- t(3)
        1457             :     !
        1458             :     !    ========[x'w'(2) = -K_zm(2)*d(xm)/dz]===rho_ds_zm(2)========== m(2)
        1459             :     !
        1460             :     !    --xm(2)---invrs_rho_ds_zt(2)---d[rho_ds_zm*K_zm*d(xm)/dz]/dz-- t(2)
        1461             :     !
        1462             :     !    ========[x'w'|_sfc]=====================rho_ds_zm(1)========== m(1) sfc
        1463             :     !
        1464             :     !    --xm(1)-------(below surface; not applicable)----------------- t(1)
        1465             :     !
        1466             :     !    where "sfc" is the level of the model surface or lower boundary.
        1467             :     !
        1468             :     !    The vertically discretized form of the turbulent advection term
        1469             :     !    (treated as an eddy diffusion term), with the explicit surface flux,
        1470             :     !    x'w'|_sfc, in place, is written out as:
        1471             :     !
        1472             :     !    - invrs_rho_ds_zt(2)
        1473             :     !      * dzt(2) * [ rho_ds_zm(2) * x'w'(2) - rho_ds_zm(1) * x'w'|_sfc ];
        1474             :     !
        1475             :     !    which can be re-written as:
        1476             :     !
        1477             :     !    + invrs_rho_ds_zt(2)
        1478             :     !      * dzt(2)
        1479             :     !        * [   rho_ds_zm(2) * K_zm(2) * dzm(2) * ( xm(3) - xm(2) )
        1480             :     !            + rho_ds_zm(1) * x'w'|_sfc ];
        1481             :     !
        1482             :     !    which can be re-written again as:
        1483             :     !
        1484             :     !    + invrs_rho_ds_zt(2)
        1485             :     !      * dzt(2)
        1486             :     !        * rho_ds_zm(2) * K_zm(2) * dzm(2) * ( xm(3) - xm(2) )
        1487             :     !    + invrs_rho_ds_zt(2)
        1488             :     !      * dzt(2)
        1489             :     !        * rho_ds_zm(1) * x'w'|_sfc.
        1490             :     !
        1491             :     !    For this equation, a Crank-Nicholson (semi-implicit) diffusion scheme
        1492             :     !    is used to solve the (1/rho_ds) * d [ rho_ds * K_zm * d(xm)/dz ] / dz
        1493             :     !    eddy-diffusion term.  The discretized implicit portion of the term is
        1494             :     !    written out as:
        1495             :     !
        1496             :     !    + (1/2) * invrs_rho_ds_zt(2)
        1497             :     !      * dzt(2)
        1498             :     !        * [ rho_ds_zm(2) * K_zm(2)
        1499             :     !            * dzm(2) * ( xm(3,<t+1>) - xm(2,<t+1>) ) ].
        1500             :     !
        1501             :     !    Note:  When the implicit term is brought over to the left-hand side,
        1502             :     !           the sign is reversed and the leading "+" in front of the term
        1503             :     !           is changed to a "-".
        1504             :     !
        1505             :     !    The discretized explicit portion of the term is written out as:
        1506             :     !
        1507             :     !    + (1/2) * invrs_rho_ds_zt(2)
        1508             :     !      * dzt(2)
        1509             :     !        * [ rho_ds_zm(2) * K_zm(2)
        1510             :     !            * dzm(2) * ( xm(3,<t>) - xm(2,<t>) ) ]
        1511             :     !    + invrs_rho_ds_zt(2)
        1512             :     !      * dzt(2)
        1513             :     !        * rho_ds_zm(1) * x'w'|_sfc.
        1514             :     !
        1515             :     !    Note:  The x'w'|_sfc portion of the term written above has been pulled
        1516             :     !           away from the rest of the explicit form written above because
        1517             :     !           the (1/2) factor due to Crank-Nicholson time_stepping does not
        1518             :     !           apply to it, as there isn't an implicit portion for x'w'|_sfc.
        1519             :     !
        1520             :     !    Timestep index (t) stands for the index of the current timestep, while
        1521             :     !    timestep index (t+1) stands for the index of the next timestep, which
        1522             :     !    is being advanced to in solving the d(xm)/dt equation.
        1523             :     !
        1524             :     ! 2) x'w' surface flux; implicit form for momentum fluxes u'w' and v'w'
        1525             :     !
        1526             :     !    The x'w' surface flux is applied to the d(xm)/dt equation through the
        1527             :     !    turbulent advection term, which is:
        1528             :     !
        1529             :     !    - (1/rho_ds) * d( rho_ds * x'w' )/dz.
        1530             :     !
        1531             :     !    At most vertical levels, a substitution can be made for x'w', such
        1532             :     !    that:
        1533             :     !
        1534             :     !    x'w' = - K_zm * d(xm)/dz.
        1535             :     !
        1536             :     !    However, the same substitution cannot be made at the surface (momentum
        1537             :     !    level 1), as x'w'|_sfc is a surface momentum flux that is found by the
        1538             :     !    following equation:
        1539             :     !
        1540             :     !    x'w'|_sfc = - [ u_star^2 / sqrt( um^2 + vm^2 ) ] * xm;
        1541             :     !
        1542             :     !    where x'w'|_sfc and xm are either u'w'|_sfc and um, respectively, or
        1543             :     !    v'w'|_sfc and vm, respectively (um and vm are located at the first
        1544             :     !    thermodynamic level above the surface, which is thermodynamic level 2),
        1545             :     !    sqrt( um^2 + vm^2 ) is the wind speed (also at thermodynamic level 2),
        1546             :     !    and u_star is defined as:
        1547             :     !
        1548             :     !    u_star = ( u'w'|_sfc^2 + v'w'|_sfc^2 )^(1/4);
        1549             :     !
        1550             :     !    and thus u_star^2 is defined as:
        1551             :     !
        1552             :     !    u_star^2 = sqrt( u'w'|_sfc^2 + v'w'|_sfc^2 ).
        1553             :     !
        1554             :     !    The value of u_star is either set to a constant value or computed
        1555             :     !    (through function diag_ustar) based on the surface wind speed, the
        1556             :     !    height above surface of the surface wind speed (as compared to the
        1557             :     !    roughness height), and the buoyancy flux at the surface.  Either way,
        1558             :     !    u_star is computed elsewhere in the model, and the values of u'w'|_sfc
        1559             :     !    and v'w'|_sfc are based on it and computed along with it.  The values
        1560             :     !    of u'w'|_sfc and v'w'|_sfc are then passed into advance_clubb_core,
        1561             :     !    and are eventually passed into advance_windm_edsclrm.  In subroutine
        1562             :     !    advance_windm_edsclrm, the value of u_star_sqd is then recomputed
        1563             :     !    based on u'w'|_sfc and v'w'|_sfc.  The value of sqrt( u_star_sqd ) is
        1564             :     !    consistent with the value of the original computation of u_star.
        1565             :     !
        1566             :     !    The equation listed above is substituted for x'w'|_sfc.  The lower
        1567             :     !    boundary condition, which in this case needs to be applied to the
        1568             :     !    d(xm)/dt equation at level 2, is discretized as follows:
        1569             :     !
        1570             :     !    --xm(3)------------------------------------------------------- t(3)
        1571             :     !
        1572             :     !    ===[x'w'(2) = -K_zm(2)*d(xm)/dz]=================rho_ds_zm(2)= m(2)
        1573             :     !
        1574             :     !    --xm(2)---invrs_rho_ds_zt(2)---d[rho_ds_zm*K_zm*d(xm)/dz]/dz-- t(2)
        1575             :     !
        1576             :     !    ===[x'w'|_sfc = -[u_star^2/sqrt(um^2+vm^2)]*xm]==rho_ds_zm(1)= m(1) sfc
        1577             :     !
        1578             :     !    --xm(1)-------(below surface; not applicable)----------------- t(1)
        1579             :     !
        1580             :     !    where "sfc" is the level of the model surface or lower boundary.
        1581             :     !
        1582             :     !    The vertically discretized form of the turbulent advection term
        1583             :     !    (treated as an eddy diffusion term), with the implicit surface momentum
        1584             :     !    flux in place, is written out as:
        1585             :     !
        1586             :     !    - invrs_rho_ds_zt(2)
        1587             :     !      * dzt(2) * [ rho_ds_zm(2) * x'w'(2) - rho_ds_zm(1) * x'w'|_sfc ];
        1588             :     !
        1589             :     !    which can be re-written as:
        1590             :     !
        1591             :     !    - invrs_rho_ds_zt(2)
        1592             :     !      * dzt(2)
        1593             :     !        * [   rho_ds_zm(2)
        1594             :     !              * { - K_zm(2) * dzm(2) * ( xm(3) - xm(2) ) }
        1595             :     !            - rho_ds_zm(1)
        1596             :     !              * { - [ u_star^2 / sqrt( um(2)^2 + vm(2)^2 ) ] * xm(2) } ];
        1597             :     !
        1598             :     !    which can be re-written as:
        1599             :     !
        1600             :     !    + invrs_rho_ds_zt(2)
        1601             :     !      * dzt(2)
        1602             :     !        * rho_ds_zm(2) * K_zm(2) * dzm(2) * ( xm(3) - xm(2) )
        1603             :     !    - invrs_rho_ds_zt(2)
        1604             :     !      * dzt(2)
        1605             :     !        * rho_ds_zm(1) * [ u_star^2 / sqrt( um(2)^2 + vm(2)^2 ) ] * xm(2).
        1606             :     !
        1607             :     !    For this equation, a Crank-Nicholson (semi-implicit) diffusion scheme
        1608             :     !    is used to solve the (1/rho_ds) * d [ rho_ds * K_zm * d(xm)/dz ] / dz
        1609             :     !    eddy-diffusion term.  The discretized implicit portion of the term is
        1610             :     !    written out as:
        1611             :     !
        1612             :     !    + (1/2) * invrs_rho_ds_zt(2)
        1613             :     !      * dzt(2)
        1614             :     !        * [ rho_ds_zm(2) * K_zm(2)
        1615             :     !            * dzm(2) * ( xm(3,<t+1>) - xm(2,<t+1>) ) ]
        1616             :     !    - invrs_rho_ds_zt(2)
        1617             :     !      * dzt(2)
        1618             :     !        * rho_ds_zm(1)
        1619             :     !        * [u_star^2/sqrt( um(2,<t>)^2 + vm(2,<t>)^2 )] * xm(2,<t+1>).
        1620             :     !
        1621             :     !    Note:  When the implicit term is brought over to the left-hand side,
        1622             :     !           the signs are reversed and the leading "+" in front of the first
        1623             :     !           part of the term is changed to a "-", while the leading "-" in
        1624             :     !           front of the second part of the term is changed to a "+".
        1625             :     !
        1626             :     !    Note:  The x'w'|_sfc portion of the term written above has been pulled
        1627             :     !           away from the rest of the implicit form written above because
        1628             :     !           the (1/2) factor due to Crank-Nicholson time_stepping does not
        1629             :     !           apply to it.  The x'w'|_sfc portion of the term is treated
        1630             :     !           completely implicitly in order to enhance numerical stability.
        1631             :     !
        1632             :     !    The discretized explicit portion of the term is written out as:
        1633             :     !
        1634             :     !    + (1/2) * invrs_rho_ds_zt(2)
        1635             :     !      * dzt(2)
        1636             :     !        * [ rho_ds_zm(2) * K_zm(2)
        1637             :     !            * dzm(2) * ( xm(3,<t>) - xm(2,<t>) ) ].
        1638             :     !
        1639             :     !    Timestep index (t) stands for the index of the current timestep, while
        1640             :     !    timestep index (t+1) stands for the index of the next timestep, which
        1641             :     !    is being advanced to in solving the d(xm)/dt equation.
        1642             :     !
        1643             :     !
        1644             :     ! The lower boundary condition for the implicit and explicit portions of the
        1645             :     ! turbulent advection term, without the x'w'|_sfc portion of the term, can
        1646             :     ! easily be invoked by using the zero-flux boundary conditions found in the
        1647             :     ! generalized diffusion function (function diffusion_zt_lhs), which is used
        1648             :     ! for many other equations in this model.  Either the generalized explicit
        1649             :     ! surface flux needs to be added onto the explicit term after the diffusion
        1650             :     ! function has been called from subroutine windm_edsclrm_rhs, or the
        1651             :     ! implicit momentum surface flux needs to be added onto the implicit term
        1652             :     ! after the diffusion function has been called from subroutine
        1653             :     ! windm_edsclrm_lhs.  However, all other equations in this model that use
        1654             :     ! zero-flux diffusion have level 1 as the level to which the lower boundary
        1655             :     ! condition needs to be applied.  Thus, an adjuster will have to be used at
        1656             :     ! level 2 to call diffusion_zt_lhs with level 1 as the input level (the last
        1657             :     ! variable being passed in during the function call).  However, the other
        1658             :     ! variables passed in (rho_ds_zm*K_zm, gr%dzt(1,:), and gr%dzm(1,:) variables) will
        1659             :     ! have to be passed in as solving for level 2.
        1660             :     !
        1661             :     ! The value of xm(1) is located below the model surface and does not effect
        1662             :     ! the rest of the model.  Since xm can be either a horizontal wind component
        1663             :     ! or a generic eddy scalar quantity, the value of xm(1) is simply set to the
        1664             :     ! value of xm(2) after the equation matrix has been solved.
        1665             :     !
        1666             :     !
        1667             :     ! Conservation Properties:
        1668             :     !
        1669             :     ! When a fixed-flux lower boundary condition is used (combined with a
        1670             :     ! zero-flux upper boundary condition), this technique of discretizing the
        1671             :     ! turbulent advection term (treated as an eddy-diffusion term) leads to
        1672             :     ! conservative differencing.  When the implicit momentum surface flux is
        1673             :     ! either zero or not used, the column totals for each column in the
        1674             :     ! left-hand side matrix (for the turbulent advection term) should be equal
        1675             :     ! to 0.  Otherwise, the column total for the second column will be equal to
        1676             :     ! rho_ds_zm(1) * x'w'|_sfc<t+1>.   When the generalized explicit surface
        1677             :     ! flux is either zero or not used, the column total for the right-hand side
        1678             :     ! vector (for the turbulent advection term) should be equal to 0.
        1679             :     ! Otherwise, the column total for the right-hand side vector (for the
        1680             :     ! turbulent advection term) will be equal to rho_ds_zm(1) * x'w'|_sfc<t>.
        1681             :     ! This ensures that the total amount of quantity xm over the entire vertical
        1682             :     ! domain is only changed by the surface flux (neglecting any forcing terms).
        1683             :     ! The total amount of change is equal to rho_ds_zm(1) * x'w'|_sfc.
        1684             :     !
        1685             :     ! To see that this conservation law is satisfied by the left-hand side
        1686             :     ! matrix, compute the turbulent advection (treated as eddy diffusion) of xm,
        1687             :     ! neglecting any implicit momentum surface flux, multiply by rho_ds_zt, and
        1688             :     ! integrate vertically.  In discretized matrix notation (where "i" stands
        1689             :     ! for the matrix column and "j" stands for the matrix row):
        1690             :     !
        1691             :     !  0 = Sum_j Sum_i
        1692             :     !       (rho_ds_zt)_i ( 1/dzt )_i
        1693             :     !       ( 0.5_core_rknd * (1/rho_ds_zt) * dzt * (rho_ds_zm*K_zm*dzm) )_ij (xm<t+1>)_j.
        1694             :     !
        1695             :     ! The left-hand side matrix,
        1696             :     ! ( 0.5_core_rknd * (1/rho_ds_zt) * dzt * (rho_ds_zm*K_zm*dzm) )_ij, is partially
        1697             :     ! written below.  The sum over i in the above equation removes (1/rho_ds_zt)
        1698             :     ! and dzt everywhere from the matrix below.  The sum over j leaves the
        1699             :     ! column totals that are desired, which are 0.
        1700             :     !
        1701             :     ! Left-hand side matrix contributions from the turbulent advection term
        1702             :     ! (treated as an eddy-diffusion term using a Crank-Nicholson timestep);
        1703             :     ! first five vertical levels:
        1704             :     !
        1705             :     !     ------------------------------------------------------------------------------->
        1706             :     !k=1 |  0             0                        0                          0
        1707             :     !    |
        1708             :     !k=2 |  0   +0.5*                  -0.5*                                  0
        1709             :     !    |        (1/rho_ds_zt(k))*      (1/rho_ds_zt(k))*
        1710             :     !    |        dzt(k)*                dzt(k)*
        1711             :     !    |        rho_ds_zm(k)*          rho_ds_zm(k)*
        1712             :     !    |        K_zm(k)*dzm(k)         K_zm(k)*dzm(k)
        1713             :     !    |
        1714             :     !k=3 |  0   -0.5*                  +0.5*                      -0.5*
        1715             :     !    |        (1/rho_ds_zt(k))*      (1/rho_ds_zt(k))*          (1/rho_ds_zt(k))*
        1716             :     !    |        dzt(k)*                dzt(k)*                    dzt(k)*
        1717             :     !    |        rho_ds_zm(k-1)*        [ rho_ds_zm(k)*            rho_ds_zm(k)*
        1718             :     !    |        K_zm(k-1)*dzm(k-1)       K_zm(k)*dzm(k)           K_zm(k)*dzm(k)
        1719             :     !    |                                +rho_ds_zm(k-1)*
        1720             :     !    |                                 K_zm(k-1)*dzm(k-1) ]
        1721             :     !    |
        1722             :     !k=4 |  0             0            -0.5*                      +0.5*
        1723             :     !    |                               (1/rho_ds_zt(k))*          (1/rho_ds_zt(k))*
        1724             :     !    |                               dzt(k)*                    dzt(k)*
        1725             :     !    |                               rho_ds_zm(k-1)*            [ rho_ds_zm(k)*
        1726             :     !    |                               K_zm(k-1)*dzm(k-1)           K_zm(k)*dzm(k)
        1727             :     !    |                                                           +rho_ds_zm(k-1)*
        1728             :     !    |                                                            K_zm(k-1)*dzm(k-1) ]
        1729             :     !    |
        1730             :     !k=5 |  0             0                        0              -0.5*
        1731             :     !    |                                                          (1/rho_ds_zt(k))*
        1732             :     !    |                                                          dzt(k)*
        1733             :     !    |                                                          rho_ds_zm(k-1)*
        1734             :     !    |                                                          K_zm(k-1)*dzm(k-1)
        1735             :     !   \ /
        1736             :     !
        1737             :     ! Note:  The superdiagonal term from level 4 and both the main diagonal and
        1738             :     !        superdiagonal terms from level 5 are not shown on this diagram.
        1739             :     !
        1740             :     ! Note:  If an implicit momentum surface flux is used, an additional term,
        1741             :     !        + (1/rho_ds_zt(2)) * dzt(2) * rho_ds_zm(1)
        1742             :     !          * [ u_star^2 / sqrt( um(2,<t>)^2 + vm(2,<t>)^2 ) ], is added to
        1743             :     !        row 2 (k=2), column 2.
        1744             :     !
        1745             :     ! To see that the above conservation law is satisfied by the right-hand side
        1746             :     ! vector, compute the turbulent advection (treated as eddy diffusion) of xm,
        1747             :     ! neglecting any generalized explicit surface flux, multiply by rho_ds_zt,
        1748             :     ! and integrate vertically.  In discretized matrix notation (where "i"
        1749             :     ! stands for the matrix column and "j" stands for the matrix row):
        1750             :     !
        1751             :     !  0 = Sum_j Sum_i (rho_ds_zt)_i ( 1/dzt )_i ( rhs_vector )_j.
        1752             :     !
        1753             :     ! The right-hand side vector, ( rhs_vector )_j, is partially written below.
        1754             :     ! The sum over i in the above equation removes (1/rho_ds_zt) and dzt
        1755             :     ! everywhere from the vector below.  The sum over j leaves the column total
        1756             :     ! that is desired, which is 0.
        1757             :     !
        1758             :     ! Right-hand side vector contributions from the turbulent advection term
        1759             :     ! (treated as an eddy-diffusion term using a Crank-Nicholson timestep);
        1760             :     ! first five vertical levels:
        1761             :     !
        1762             :     !     --------------------------------------------
        1763             :     !k=1 |                      0                     |
        1764             :     !    |                                            |
        1765             :     !    |                                            |
        1766             :     !k=2 | +0.5*(1/rho_ds_zt(k))*                     |
        1767             :     !    |      dzt(k)*                               |
        1768             :     !    |       [ rho_ds_zm(k)*K_zm(k)*              |
        1769             :     !    |         dzm(k)*(xm(k+1,<t>)-xm(k,<t>)) ]   |
        1770             :     !    |                                            |
        1771             :     !k=3 | +0.5*(1/rho_ds_zt(k))*                     |
        1772             :     !    |      dzt(k)*                               |
        1773             :     !    |       [ rho_ds_zm(k)*K_zm(k)*              |
        1774             :     !    |         dzm(k)*(xm(k+1,<t>)-xm(k,<t>))     |
        1775             :     !    |        -rho_ds_zm(k-1)*K_zm(k-1)*          |
        1776             :     !    |         dzm(k-1)*(xm(k,<t>)-xm(k-1,<t>)) ] |
        1777             :     !    |                                            |
        1778             :     !k=4 | +0.5*(1/rho_ds_zt(k))*                     |
        1779             :     !    |      dzt(k)*                               |
        1780             :     !    |       [ rho_ds_zm(k)*K_zm(k)*              |
        1781             :     !    |         dzm(k)*(xm(k+1,<t>)-xm(k,<t>))     |
        1782             :     !    |        -rho_ds_zm(k-1)*K_zm(k-1)*          |
        1783             :     !    |         dzm(k-1)*(xm(k,<t>)-xm(k-1,<t>)) ] |
        1784             :     !    |                                            |
        1785             :     !k=5 | +0.5*(1/rho_ds_zt(k))*                     |
        1786             :     !    |      dzt(k)*                               |
        1787             :     !    |       [ rho_ds_zm(k)*K_zm(k)*              |
        1788             :     !    |         dzm(k)*(xm(k+1,<t>)-xm(k,<t>))     |
        1789             :     !    |        -rho_ds_zm(k-1)*K_zm(k-1)*          |
        1790             :     !    |         dzm(k-1)*(xm(k,<t>)-xm(k-1,<t>)) ] |
        1791             :     !   \ /                                          \ /
        1792             :     !
        1793             :     ! Note:  If a generalized explicit surface flux is used, an additional term,
        1794             :     !        + (1/rho_ds_zt(2)) * dzt(2) * rho_ds_zm(1) * x'w'|_sfc, is added to
        1795             :     !        row 2 (k=2).
        1796             :     !
        1797             :     ! Note:  Only the contributions by the turbulent advection term are shown
        1798             :     !        for both the left-hand side matrix and the right-hand side vector.
        1799             :     !        There are more terms in the equation, and thus more factors to be
        1800             :     !        added to both the left-hand side matrix (such as time tendency and
        1801             :     !        mean advection) and the right-hand side vector (such as xm
        1802             :     !        forcings).  The left-hand side matrix is set-up so that a singular
        1803             :     !        matrix is not encountered.
        1804             :
        1805             :     ! References:
        1806             :     ! Eqn. 8 & 9 on p. 3545 of
        1807             :     ! ``A PDF-Based Model for Boundary Layer Clouds. Part I:
        1808             :     !   Method and Model Description'' Golaz, et al. (2002)
        1809             :     ! JAS, Vol. 59, pp. 3540--3551.
        1810             :     !-----------------------------------------------------------------------
        1811             :
        1812             :     use grid_class, only: &
        1813             :         grid ! Type
        1814             :
        1815             :     use matrix_solver_wrapper, only:  &
        1816             :         tridiag_solve ! Procedure(s)
        1817             :
        1818             :     use stats_variables, only: &
        1819             :         stats_metadata_type
        1820             :
        1821             :     use stats_type_utilities, only:  &
        1822             :         stat_update_var_pt  ! Subroutine
        1823             :
        1824             :     use clubb_precision, only: &
        1825             :         core_rknd ! Variable(s)
        1826             :
        1827             :     use stats_type, only: stats ! Type
        1828             :
        1829             :     implicit none
        1830             :
        1831             :     ! Constant parameters
        1832             :
        1833             :     integer, parameter :: &
        1834             :       kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
        1835             :       k_tdiag   = 2, & ! Thermodynamic main diagonal index.
        1836             :       km1_tdiag = 3    ! Thermodynamic subdiagonal index.
        1837             :
        1838             :     ! ------------------------ Input Variables ------------------------
        1839             :     integer, intent(in) :: &
        1840             :       nz, &
        1841             :       ngrdcol
        1842             :
        1843             :     type (grid), target, intent(in) :: &
        1844             :       gr
        1845             :
        1846             :     integer, intent(in) :: &
        1847             :       nrhs ! Number of right-hand side (explicit) vectors & Number of solution vectors.
        1848             :
        1849             :     integer, intent(in) :: &
        1850             :       ixm_matrix_condt_num  ! Stats index of the condition numbers
        1851             :
        1852             :     integer, intent(in) :: &
        1853             :       tridiag_solve_method  ! Specifier for method to solve tridiagonal systems
        1854             :
        1855             :     type (stats_metadata_type), intent(in) :: &
        1856             :       stats_metadata
        1857             :
        1858             :     ! ------------------------ Inout variables ------------------------
        1859             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
        1860             :       stats_sfc
        1861             :
        1862             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(inout) :: &
        1863             :       lhs    ! Implicit contributions to um, vm, and eddy scalars  [units vary]
        1864             :
        1865             :     real( kind = core_rknd ), dimension(ngrdcol,nz,nrhs), intent(inout) :: &
        1866             :       rhs    ! Right-hand side (explicit) contributions.
        1867             :
        1868             :     ! ------------------------ Output variables ------------------------
        1869             :     real( kind = core_rknd ), dimension(ngrdcol,nz,nrhs), intent(out) :: &
        1870             :       solution ! Solution to the system of equations    [units vary]
        1871             :
        1872             :     ! ------------------------ Local variables ------------------------
        1873             :     real( kind = core_rknd ), dimension(ngrdcol) :: &
        1874    17870112 :       rcond ! Estimate of the reciprocal of the condition number on the LHS matrix
        1875             :
        1876             :     integer :: i
        1877             :
        1878             :     ! ------------------------ Begin Code ------------------------
        1879             :
        1880             :     ! Solve tridiagonal system for xm.
        1881     8935056 :     if ( stats_metadata%l_stats_samp .and. ixm_matrix_condt_num > 0 ) then
        1882             :
        1883             :       call tridiag_solve( "windm_edsclrm", tridiag_solve_method,  & ! Intent(in)
        1884             :                           ngrdcol, nz, nrhs,                      & ! Intent(in)
        1885             :                           lhs, rhs,                               & ! Intent(inout)
        1886           0 :                           solution, rcond )                         ! Intent(out)
        1887             :
        1888             :       ! Est. of the condition number of the variance LHS matrix
        1889           0 :       do i = 1, ngrdcol
        1890           0 :         call stat_update_var_pt( ixm_matrix_condt_num, 1, 1.0_core_rknd/rcond(i), &  ! intent(in)
        1891           0 :                                  stats_sfc(i) )                                      ! intent(inout)
        1892             :       end do
        1893             :     else
        1894             :
        1895             :       call tridiag_solve( "windm_edsclrm", tridiag_solve_method,  & ! Intent(in)
        1896             :                           ngrdcol, nz, nrhs,                      & ! Intent(in)
        1897             :                           lhs, rhs,                               & ! Intent(inout)
        1898     8935056 :                           solution )                                ! Intent(out)
        1899             :     end if
        1900             :
        1901     8935056 :     return
        1902             :   end subroutine windm_edsclrm_solve
        1903             :
        1904             :   !=============================================================================
        1905           0 :   subroutine windm_edsclrm_implicit_stats( nz, solve_type, xm, &
        1906           0 :                                            invrs_dzt, &
        1907           0 :                                            lhs_diff, lhs_ma_zt, &
        1908           0 :                                            invrs_rho_ds_zt, u_star_sqd,&
        1909           0 :                                            rho_ds_zm, wind_speed, &
        1910             :                                            l_imp_sfc_momentum_flux, &
        1911             :                                            stats_metadata, &
        1912             :                                            stats_zt )
        1913             :
        1914             :     ! Description:
        1915             :     ! Compute implicit contributions to um and vm
        1916             :
        1917             :     ! References:
        1918             :     ! None
        1919             :     !-----------------------------------------------------------------------
        1920             :
        1921             :     use constants_clubb, only: &
        1922             :         zero
        1923             :
        1924             :     use stats_type_utilities, only:  &
        1925             :         stat_end_update_pt,  & ! Subroutines
        1926             :         stat_update_var_pt
        1927             :
        1928             :     use clubb_precision, only:  &
        1929             :         core_rknd
        1930             :
        1931             :     use stats_type, only: &
        1932             :         stats ! Type
        1933             :
        1934             :     use stats_variables, only: &
        1935             :         stats_metadata_type
        1936             :
        1937             :     implicit none
        1938             :
        1939             :     !---------------------- Input Variables ----------------------
        1940             :     integer, intent(in) :: &
        1941             :       nz
        1942             :
        1943             :     integer, intent(in) :: &
        1944             :       solve_type     ! Desc. of what is being solved for
        1945             :
        1946             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
        1947             :       xm,         & !  Computed value um or vm at <t+1>    [m/s]
        1948             :       invrs_dzt     ! The inverse spacing between momentum grid levels;
        1949             :                     ! centered over thermodynamic grid levels.
        1950             :
        1951             :     real( kind = core_rknd ), dimension(3,nz), intent(in) :: &
        1952             :       lhs_diff, & ! LHS diffustion terms
        1953             :       lhs_ma_zt   ! LHS mean advection terms
        1954             :
        1955             :     real( kind = core_rknd ), dimension(nz), intent(in) ::  &
        1956             :       wind_speed,      & ! wind speed; sqrt(u^2 + v^2)              [m/s]
        1957             :       rho_ds_zm,       & ! Dry, static density on momentum levels      [kg/m^3]
        1958             :       invrs_rho_ds_zt    ! Inv. dry, static density at thermo. levels  [m^3/kg]
        1959             :
        1960             :     real( kind = core_rknd ), intent(in) :: &
        1961             :       u_star_sqd   ! Surface friction velocity, u_star, squared      [m/s]
        1962             :
        1963             :     logical, intent(in) :: &
        1964             :       l_imp_sfc_momentum_flux  ! Flag for implicit momentum surface fluxes.
        1965             :
        1966             :     type (stats_metadata_type), intent(in) :: &
        1967             :       stats_metadata
        1968             :
        1969             :     !---------------------- InOut variables ----------------------
        1970             :     type (stats), target, intent(inout) :: &
        1971             :       stats_zt
        1972             :
        1973             :     !---------------------- Local variables ----------------------
        1974             :     integer :: k, kp1, km1 ! Array indices
        1975             :
        1976             :     real( kind = core_rknd ), dimension(nz) :: &
        1977           0 :       imp_sfc_flux
        1978             :
        1979             :     ! Budget indices
        1980             :     integer :: ixm_ma, ixm_ta
        1981             :
        1982             :     !---------------------- Begin Code ----------------------
        1983             :
        1984           0 :     select case ( solve_type )
        1985             :     case ( windm_edsclrm_um )
        1986           0 :       ixm_ma = stats_metadata%ium_ma
        1987           0 :       ixm_ta = stats_metadata%ium_ta
        1988             :
        1989             :     case ( windm_edsclrm_vm )
        1990           0 :       ixm_ma = stats_metadata%ivm_ma
        1991           0 :       ixm_ta = stats_metadata%ivm_ta
        1992             :
        1993             :     case default
        1994           0 :       ixm_ma = 0
        1995           0 :       ixm_ta = 0
        1996             :
        1997             :     end select
        1998             :
        1999           0 :     imp_sfc_flux(:) = zero
        2000             :
        2001           0 :     if ( l_imp_sfc_momentum_flux ) then
        2002             :
        2003             :       ! Statistics:  implicit contributions for um or vm.
        2004             :
        2005             :       ! xm term ta is modified at level 2 to include the effects of the
        2006             :       ! surface flux.  In this case, this effects the implicit portion of
        2007             :       ! the term, which handles the main diagonal for the turbulent advection term
        2008           0 :       if ( stats_metadata%ium_ta + stats_metadata%ivm_ta > 0 ) then
        2009           0 :         imp_sfc_flux(2) =  - invrs_rho_ds_zt(2) * invrs_dzt(2) &
        2010           0 :                              * rho_ds_zm(1) * ( u_star_sqd / wind_speed(2) )
        2011             :       endif
        2012             :
        2013             :     endif ! l_imp_sfc_momentum_flux
        2014             :
        2015             :     ! Finalize implicit contributions for xm
        2016             :
        2017           0 :     do k = 2, nz-1, 1
        2018             :
        2019           0 :       km1 = max( k-1, 2 )
        2020           0 :       kp1 = min( k+1, nz )
        2021             :
        2022             :       ! xm mean advection
        2023             :       ! xm term ma is completely implicit; call stat_update_var_pt.
        2024             :       call stat_update_var_pt( ixm_ma, k, & ! intent(in)
        2025           0 :              - lhs_ma_zt(3,k) * xm(km1)   &
        2026             :              - lhs_ma_zt(2,k) * xm(k)     &
        2027           0 :              - lhs_ma_zt(1,k) * xm(kp1),  & ! intent(in)
        2028           0 :               stats_zt )                    ! intent(inout)
        2029             :
        2030             :       ! xm turbulent transport (implicit component)
        2031             :       ! xm term ta has both implicit and explicit components;
        2032             :       ! call stat_end_update_pt.
        2033             :       call stat_end_update_pt( ixm_ta, k,               & ! intent(in)
        2034           0 :              - 0.5_core_rknd * lhs_diff(3,k) * xm(km1)  &
        2035             :              + ( - 0.5_core_rknd * lhs_diff(2,k)        &
        2036             :                  + imp_sfc_flux(k) ) * xm(k)               &
        2037             :              - 0.5_core_rknd * lhs_diff(1,k) * xm(kp1), & ! intent(in)
        2038           0 :            stats_zt )                                     ! intent(inout)
        2039             :
        2040             :     enddo
        2041             :
        2042             :
        2043             :     ! Upper boundary conditions
        2044           0 :     k   = nz
        2045           0 :     km1 = max( k-1, 1 )
        2046             :
        2047             :     ! xm mean advection
        2048             :     ! xm term ma is completely implicit; call stat_update_var_pt.
        2049             :     call stat_update_var_pt( ixm_ma, k, & ! intent(in)
        2050           0 :            - lhs_ma_zt(3,k) * xm(km1)   &
        2051             :            - lhs_ma_zt(2,k) * xm(k),    & ! intent(in)
        2052           0 :             stats_zt )                    ! intent(inout)
        2053             :
        2054             :     ! xm turbulent transport (implicit component)
        2055             :     ! xm term ta has both implicit and explicit components;
        2056             :     ! call stat_end_update_pt.
        2057             :     call stat_end_update_pt( ixm_ta, k,               & ! intent(in)
        2058           0 :            - 0.5_core_rknd * lhs_diff(3,k) * xm(km1)  &
        2059             :            + ( - 0.5_core_rknd * lhs_diff(2,k)        &
        2060             :                + imp_sfc_flux(k) ) * xm(k),              & ! intent(in)
        2061           0 :            stats_zt )                                   ! intent(inout)
        2062             :
        2063             :
        2064           0 :     return
        2065             :   end subroutine windm_edsclrm_implicit_stats
        2066             :
        2067             :   !=============================================================================
        2068           0 :   subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
        2069           0 :                                fcor, perp_wind_m, perp_wind_g, &
        2070           0 :                                xm_forcing, l_implemented, &
        2071             :                                stats_metadata, &
        2072           0 :                                stats_zt, &
        2073           0 :                                xm_tndcy )
        2074             :
        2075             :     ! Description:
        2076             :     ! Computes the explicit tendency for the um and vm wind components.
        2077             :     !
        2078             :     ! The only explicit tendency that is involved in the d(um)/dt or d(vm)/dt
        2079             :     ! equations is the Coriolis tendency.
        2080             :     !
        2081             :     ! The d(um)/dt equation contains the term:
        2082             :     !
        2083             :     ! - f * ( v_g - vm );
        2084             :     !
        2085             :     ! where f is the Coriolis parameter and v_g is the v component of the
        2086             :     ! geostrophic wind.
        2087             :     !
        2088             :     ! Likewise, the d(vm)/dt equation contains the term:
        2089             :     !
        2090             :     ! + f * ( u_g - um );
        2091             :     !
        2092             :     ! where u_g is the u component of the geostrophic wind.
        2093             :     !
        2094             :     ! This term is treated completely explicitly.  The values of um, vm, u_g,
        2095             :     ! and v_g are all found on the thermodynamic levels.
        2096             :     !
        2097             :     ! Wind forcing from the GCSS cases is also added here.
        2098             :     !
        2099             :     ! References:
        2100             :     !-----------------------------------------------------------------------
        2101             :
        2102             :     use grid_class, only: &
        2103             :         grid ! Type
        2104             :
        2105             :     use stats_type_utilities, only: &
        2106             :         stat_update_var
        2107             :
        2108             :     use stats_variables, only: &
        2109             :         stats_metadata_type
        2110             :
        2111             :     use clubb_precision, only: &
        2112             :         core_rknd ! Variable(s)
        2113             :
        2114             :     use stats_type, only: stats ! Type
        2115             :
        2116             :     implicit none
        2117             :
        2118             :     ! -------------------------- Input Variables --------------------------
        2119             :     integer, intent(in) :: &
        2120             :       nz, &
        2121             :       ngrdcol
        2122             :
        2123             :     integer, intent(in) ::  &
        2124             :       solve_type      ! Description of what is being solved for
        2125             :
        2126             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
        2127             :       fcor            ! Coriolis parameter     [s^-1]
        2128             :
        2129             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        2130             :       perp_wind_m,  & ! Perpendicular component of the mean wind (e.g. v, for the u-eqn) [m/s]
        2131             :       perp_wind_g,  & ! Perpendicular component of the geostropic wind (e.g. vg)         [m/s]
        2132             :       xm_forcing      ! Prescribed wind forcing                                          [m/s/s]
        2133             :
        2134             :     logical, intent(in) :: &
        2135             :       l_implemented   ! Flag for CLUBB being implemented in a larger model.
        2136             :
        2137             :     type (stats_metadata_type), intent(in) :: &
        2138             :       stats_metadata
        2139             :
        2140             :     ! -------------------------- InOut Variables --------------------------
        2141             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
        2142             :       stats_zt
        2143             :
        2144             :     ! -------------------------- Output Variables --------------------------
        2145             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
        2146             :       xm_tndcy        ! xm tendency            [m/s^2]
        2147             :
        2148             :     ! -------------------------- Local Variables --------------------------
        2149             :     integer :: &
        2150             :       ixm_gf, &
        2151             :       ixm_cf, &
        2152             :       ixm_f
        2153             :
        2154             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
        2155           0 :       xm_gf, &
        2156           0 :       xm_cf
        2157             :
        2158             :     integer :: i, k
        2159             :
        2160             :     ! -------------------------- Begin Code --------------------------
        2161             :
        2162             :     !$acc enter data create( xm_gf, xm_cf )
        2163             :
        2164             :     ! Initialize variables to 0.
        2165           0 :     xm_gf = 0.0_core_rknd
        2166           0 :     xm_cf = 0.0_core_rknd
        2167           0 :     xm_tndcy = 0.0_core_rknd
        2168             :
        2169           0 :     if ( .not. l_implemented ) then
        2170             :       ! Only compute the Coriolis term if the model is running on it's own,
        2171             :       ! and is not part of a larger, host model.
        2172             :
        2173           0 :       select case ( solve_type )
        2174             :
        2175             :       case ( windm_edsclrm_um )
        2176             :
        2177           0 :         ixm_gf = stats_metadata%ium_gf
        2178           0 :         ixm_cf = stats_metadata%ium_cf
        2179           0 :         ixm_f  = stats_metadata%ium_f
        2180             :
        2181             :         !$acc parallel loop gang vector collapse(2) default(present)
        2182           0 :         do k = 2, nz
        2183           0 :           do i = 1, ngrdcol
        2184           0 :             xm_gf(i,k) = - fcor(i) * perp_wind_g(i,k)
        2185             :           end do
        2186             :         end do
        2187             :         !$acc end parallel loop
        2188             :
        2189             :         !$acc parallel loop gang vector collapse(2) default(present)
        2190           0 :         do k = 2, nz
        2191           0 :           do i = 1, ngrdcol
        2192           0 :             xm_cf(i,k) = fcor(i) * perp_wind_m(i,k)
        2193             :           end do
        2194             :         end do
        2195             :         !$acc end parallel loop
        2196             :
        2197             :       case ( windm_edsclrm_vm )
        2198             :
        2199           0 :         ixm_gf = stats_metadata%ivm_gf
        2200           0 :         ixm_cf = stats_metadata%ivm_cf
        2201           0 :         ixm_f  = stats_metadata%ivm_f
        2202             :
        2203             :         !$acc parallel loop gang vector collapse(2) default(present)
        2204           0 :         do k = 2, nz
        2205           0 :           do i = 1, ngrdcol
        2206           0 :             xm_gf(i,k) = fcor(i) * perp_wind_g(i,k)
        2207             :           end do
        2208             :         end do
        2209             :         !$acc end parallel loop
        2210             :
        2211             :         !$acc parallel loop gang vector collapse(2) default(present)
        2212           0 :         do k = 2, nz
        2213           0 :           do i = 1, ngrdcol
        2214           0 :             xm_cf(i,k) = -fcor(i) * perp_wind_m(i,k)
        2215             :           end do
        2216             :         end do
        2217             :         !$acc end parallel loop
        2218             :
        2219             :       case default
        2220             :
        2221           0 :         ixm_gf = 0
        2222           0 :         ixm_cf = 0
        2223           0 :         ixm_f = 0
        2224             :
        2225             :         !$acc parallel loop gang vector collapse(2) default(present)
        2226           0 :         do k = 2, nz
        2227           0 :           do i = 1, ngrdcol
        2228           0 :             xm_gf(i,k) = 0._core_rknd
        2229           0 :             xm_cf(i,k) = 0._core_rknd
        2230             :           end do
        2231             :         end do
        2232             :         !$acc end parallel loop
        2233             :
        2234             :       end select
        2235             :
        2236             :       !$acc parallel loop gang vector collapse(2) default(present)
        2237           0 :       do k = 2, nz
        2238           0 :         do i = 1, ngrdcol
        2239           0 :           xm_tndcy(i,k) = xm_gf(i,k) + xm_cf(i,k) + xm_forcing(i,k)
        2240             :         end do
        2241             :       end do
        2242             :       !$acc end parallel loop
        2243             :
        2244           0 :       if ( stats_metadata%l_stats_samp ) then
        2245             :
        2246             :         !$acc update host( xm_gf, xm_cf, xm_forcing )
        2247             :
        2248           0 :         do i = 1, ngrdcol
        2249             :           ! xm term gf is completely explicit; call stat_update_var.
        2250           0 :           call stat_update_var( ixm_gf, xm_gf(i,:), & ! intent(in)
        2251           0 :                                 stats_zt(i) )         ! intent(inout)
        2252             :
        2253             :           ! xm term cf is completely explicit; call stat_update_var.
        2254             :           call stat_update_var( ixm_cf, xm_cf(i,:), & ! intent(in)
        2255           0 :                                 stats_zt(i) )         ! intent(inout)
        2256             :
        2257             :           ! xm term F
        2258             :           call stat_update_var( ixm_f, xm_forcing(i,:), & ! intent(in)
        2259           0 :                                 stats_zt(i) )             ! intent(inout)
        2260             :         end do
        2261             :       endif
        2262             :
        2263             :     else   ! implemented in a host model.
        2264             :
        2265             :       !$acc parallel loop gang vector collapse(2) default(present)
        2266           0 :       do k = 2, nz
        2267           0 :         do i = 1, ngrdcol
        2268           0 :           xm_tndcy(i,k) = 0.0_core_rknd
        2269             :         end do
        2270             :       end do
        2271             :       !$acc end parallel loop
        2272             :
        2273             :     endif
        2274             :
        2275             :     !$acc exit data delete( xm_gf, xm_cf )
        2276             :
        2277           0 :     return
        2278             :
        2279             :   end subroutine compute_uv_tndcy
        2280             :
        2281             :   !======================================================================================
        2282     8935056 :   subroutine windm_edsclrm_lhs( nz, ngrdcol, gr, dt, &
        2283     8935056 :                                 lhs_ma_zt, lhs_diff, &
        2284     8935056 :                                 wind_speed, u_star_sqd,  &
        2285     8935056 :                                 rho_ds_zm, invrs_rho_ds_zt,  &
        2286             :                                 l_implemented, l_imp_sfc_momentum_flux,  &
        2287     8935056 :                                 lhs )
        2288             :     ! Description:
        2289             :     ! Calculate the implicit portion of the horizontal wind or eddy-scalar
        2290             :     ! time-tendency equation.  See the description in subroutine
        2291             :     ! windm_edsclrm_solve for more details.
        2292             :     !
        2293             :     ! Notes:
        2294             :     !   Lower Boundary:
        2295             :     !       The lower boundary condition is a fixed-flux boundary condition, which
        2296             :     !       gets added into the time-tendency equation at level 2.
        2297             :     !       The value of xm(1) is located below the model surface and does not effect
        2298             :     !       the rest of the model.  Since xm can be either a horizontal wind component
        2299             :     !       or a generic eddy scalar quantity, the value of xm(1) is simply set to the
        2300             :     !       value of xm(2) after the equation matrix has been solved.
        2301             :     !
        2302             :     !   --- THIS SUBROUTINE HAS BEEN OPTIMIZED ---
        2303             :     !   Simple changes to this procedure may adversely affect computational speed
        2304             :     !       - Gunther Huebler, Aug. 2018, clubb:ticket:834
        2305             :     !----------------------------------------------------------------------------------
        2306             :
        2307             :     use grid_class, only:  &
        2308             :         grid   ! Type
        2309             :
        2310             :     use clubb_precision, only:  &
        2311             :         core_rknd ! Variable(s)
        2312             :
        2313             :     implicit none
        2314             :
        2315             :     ! ----------------------- Input Variables -----------------------
        2316             :     integer, intent(in) :: &
        2317             :       nz, &
        2318             :       ngrdcol
        2319             :
        2320             :     type (grid), target, intent(in) :: gr
        2321             :
        2322             :     real( kind = core_rknd ), intent(in) :: &
        2323             :       dt                 ! Model timestep                             [s]
        2324             :
        2325             :     real( kind = core_rknd ), intent(in), dimension(3,ngrdcol,nz) :: &
        2326             :       lhs_diff, & ! LHS diffustion terms
        2327             :       lhs_ma_zt   ! LHS mean advection terms
        2328             :
        2329             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        2330             :       wind_speed,      & ! wind speed; sqrt( u^2 + v^2 )              [m/s]
        2331             :       rho_ds_zm,       & ! Dry, static density on momentum levels     [kg/m^3]
        2332             :       invrs_rho_ds_zt    ! Inv. dry, static density at thermo. levels [m^3/kg]
        2333             :
        2334             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
        2335             :       u_star_sqd    ! Surface friction velocity, u_*, squared  [m/s]
        2336             :
        2337             :     logical, intent(in) ::  &
        2338             :       l_implemented, & ! Flag for CLUBB being implemented in a larger model.
        2339             :       l_imp_sfc_momentum_flux  ! Flag for implicit momentum surface fluxes.
        2340             :
        2341             :     ! ----------------------- Output Variable -----------------------
        2342             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(out) :: &
        2343             :       lhs           ! Implicit contributions to xm (tridiagonal matrix)
        2344             :
        2345             :     ! ----------------------- Local Variables -----------------------
        2346             :     real( kind = core_rknd ) :: &
        2347             :         invrs_dt    ! Inverse of dt, 1/dt, used for computational efficiency
        2348             :
        2349             :     integer :: k, i ! Loop variable
        2350             :
        2351             :     ! ----------------------- Begin Code -----------------------
        2352             :
        2353             :     ! Calculate coefs of eddy diffusivity and inverse of dt
        2354     8935056 :     invrs_dt = 1.0_core_rknd / dt
        2355             :
        2356             :     ! Set lower boundary, see notes
        2357             :     !$acc parallel loop gang vector default(present)
        2358   149194656 :     do i = 1, ngrdcol
        2359   140259600 :       lhs(1,i,1) = 0.0_core_rknd
        2360   140259600 :       lhs(2,i,1) = 1.0_core_rknd
        2361   149194656 :       lhs(3,i,1) = 0.0_core_rknd
        2362             :     end do
        2363             :     !$acc end parallel loop
        2364             :
        2365             :     ! Add terms to lhs
        2366             :     !$acc parallel loop gang vector collapse(2) default(present)
        2367   759479760 :     do k = 2, nz
        2368 12541286160 :       do i = 1, ngrdcol
        2369 11781806400 :         lhs(1,i,k) = 0.5_core_rknd * lhs_diff(1,i,k)
        2370             :
        2371 11781806400 :         lhs(2,i,k) = 0.5_core_rknd * lhs_diff(2,i,k)
        2372             :
        2373             :         ! LHS time tendency.
        2374 11781806400 :         lhs(2,i,k) = lhs(2,i,k) + invrs_dt
        2375             :
        2376 12532351104 :         lhs(3,i,k) = 0.5_core_rknd * lhs_diff(3,i,k)
        2377             :       end do
        2378             :     end do
        2379             :     !$acc end parallel loop
        2380             :
        2381             :     ! LHS mean advection term.
        2382     8935056 :     if ( .not. l_implemented ) then
        2383             :       !$acc parallel loop gang vector collapse(2) default(present)
        2384           0 :       do k = 2, nz-1
        2385           0 :         do i = 1, ngrdcol
        2386           0 :           lhs(1:3,i,k) = lhs(1:3,i,k) + lhs_ma_zt(:,i,k)
        2387             :         end do
        2388             :       end do
        2389             :       !$acc end parallel loop
        2390             :
        2391             :     endif
        2392             :
        2393     8935056 :     if ( l_imp_sfc_momentum_flux ) then
        2394             :
        2395             :       ! LHS momentum surface flux.
        2396             :       !$acc parallel loop gang vector default(present)
        2397           0 :       do i = 1, ngrdcol
        2398           0 :         lhs(2,i,2) = lhs(2,i,2) + invrs_rho_ds_zt(i,2) * gr%invrs_dzt(i,2) &
        2399           0 :                                   * rho_ds_zm(i,1) * ( u_star_sqd(i) / wind_speed(i,2) )
        2400             :       end do
        2401             :       !$acc end parallel loop
        2402             :     end if ! l_imp_sfc_momentum_flux
        2403             :
        2404     8935056 :     return
        2405             :
        2406             :   end subroutine windm_edsclrm_lhs
        2407             :
        2408             :   !=============================================================================
        2409   116155728 :   subroutine windm_edsclrm_rhs( nz, ngrdcol, gr, solve_type, dt, &
        2410   116155728 :                                 lhs_diff, xm, xm_tndcy,  &
        2411   116155728 :                                 rho_ds_zm, invrs_rho_ds_zt,  &
        2412   116155728 :                                 l_imp_sfc_momentum_flux, xpwp_sfc, &
        2413             :                                 stats_metadata, &
        2414   116155728 :                                 stats_zt, &
        2415   116155728 :                                 rhs )
        2416             :     ! Description:
        2417             :     !   Calculate the explicit portion of the horizontal wind or eddy-scalar
        2418             :     !   time-tendency equation.  See the description in subroutine
        2419             :     !   windm_edsclrm_solve for more details.
        2420             :     !
        2421             :     ! References:
        2422             :     !   None
        2423             :     !
        2424             :     ! Notes:
        2425             :     !   The lower boundary condition needs to be applied here at level 2.
        2426             :     !   The lower boundary condition is a "fixed flux" boundary condition.
        2427             :     !   The coding is the same as for a zero-flux boundary condition, but with
        2428             :     !   an extra term added on the right-hand side at the boundary level.  For
        2429             :     !   the rest of the model code, a zero-flux boundary condition is applied
        2430             :     !   at level 1, and thus subroutine diffusion_zt_lhs is set-up to do that.
        2431             :     !   In order to apply the same boundary condition code here at level 2, an
        2432             :     !   adjuster needs to be used to tell diffusion_zt_lhs to use the code at
        2433             :     !   level 2 that it normally uses at level 1.
        2434             :     !
        2435             :     !   --- THIS SUBROUTINE HAS BEEN OPTIMIZED ---
        2436             :     !   Simple changes to this procedure may adversely affect computational speed
        2437             :     !       - Gunther Huebler, Aug. 2018, clubb:ticket:834
        2438             :     !----------------------------------------------------------------------------------
        2439             :
        2440             :     use clubb_precision, only:  &
        2441             :         core_rknd ! Variable(s)
        2442             :
        2443             :     use diffusion, only:  &
        2444             :         diffusion_zt_lhs    ! Procedure(s)
        2445             :
        2446             :     use stats_variables, only: &
        2447             :         stats_metadata_type
        2448             :
        2449             :     use stats_type_utilities, only: &
        2450             :         stat_begin_update_pt,  & ! Procedure(s)
        2451             :         stat_modify_pt
        2452             :
        2453             :     use grid_class, only:  &
        2454             :         grid ! Type
        2455             :
        2456             :     use stats_type, only: stats ! Type
        2457             :
        2458             :     implicit none
        2459             :
        2460             :     !------------------- Input Variables -------------------
        2461             :     integer, intent(in) :: &
        2462             :       nz, &
        2463             :       ngrdcol
        2464             :
        2465             :     type (grid), target, intent(in) :: &
        2466             :       gr
        2467             :
        2468             :     integer, intent(in) :: &
        2469             :       solve_type ! Description of what is being solved for
        2470             :
        2471             :     real( kind = core_rknd ), intent(in) :: &
        2472             :       dt                 ! Model timestep                             [s]
        2473             :
        2474             :     real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(in) :: &
        2475             :       lhs_diff   ! LHS diffustion terms
        2476             :
        2477             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        2478             :       xm,              & ! Eddy-scalar variable, xm (thermo. levels)  [units vary]
        2479             :       xm_tndcy,        & ! The explicit time-tendency acting on xm    [units vary]
        2480             :       rho_ds_zm,       & ! Dry, static density on momentum levels     [kg/m^3]
        2481             :       invrs_rho_ds_zt    ! Inv. dry, static density at thermo. levels [m^3/kg]
        2482             :
        2483             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
        2484             :       xpwp_sfc     ! x'w' at the surface                              [units vary]
        2485             :
        2486             :     logical, intent(in) :: &
        2487             :       l_imp_sfc_momentum_flux  ! Flag for implicit momentum surface fluxes.
        2488             :
        2489             :     type (stats_metadata_type), intent(in) :: &
        2490             :       stats_metadata
        2491             :
        2492             :     !------------------- Inout Variable -------------------
        2493             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
        2494             :       stats_zt
        2495             :
        2496             :     !------------------- Output Variable -------------------
        2497             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        2498             :       rhs          ! Right-hand side (explicit) contributions.
        2499             :
        2500             :     !------------------- Local Variables -------------------
        2501             :     integer :: i, k    ! Loop variable
        2502             :
        2503             :     real( kind = core_rknd ) :: invrs_dt
        2504             :
        2505             :     integer :: ixm_ta
        2506             :
        2507             :     !------------------- Begin Code -------------------
        2508             :
        2509   116155728 :     invrs_dt = 1.0_core_rknd / dt   ! Precalculate 1.0/dt to avoid redoing the divide
        2510             :
        2511   116155728 :     select case ( solve_type )
        2512             :       case ( windm_edsclrm_um )
        2513           0 :         ixm_ta = stats_metadata%ium_ta
        2514             :       case ( windm_edsclrm_vm )
        2515           0 :         ixm_ta = stats_metadata%ivm_ta
        2516             :       case default  ! Eddy scalars
        2517   116155728 :         ixm_ta = 0
        2518             :     end select
        2519             :
        2520             :     ! For purposes of the matrix equation, rhs(1) is simply set to 0.
        2521             :     !$acc parallel loop gang vector default(present)
        2522  1939530528 :     do i = 1, ngrdcol
        2523  1939530528 :       rhs(i,1) = 0.0_core_rknd
        2524             :     end do
        2525             :     !$acc end parallel loop
        2526             :
        2527             :     ! Lower boundary calculation
        2528             :     !$acc parallel loop gang vector default(present)
        2529  1939530528 :     do i = 1, ngrdcol
        2530  3646749600 :       rhs(i,2) = 0.5_core_rknd  &
        2531             :                  * ( - lhs_diff(2,i,2) * xm(i,2)        &
        2532  1823374800 :                      - lhs_diff(1,i,2) * xm(i,3) )      &
        2533             :                  + xm_tndcy(i,2)                        & ! RHS forcings
        2534  3762905328 :                  + invrs_dt * xm(i,2)                     ! RHS time tendency
        2535             :     end do
        2536             :     !$acc end parallel loop
        2537             :
        2538             :     ! Non-boundary rhs calculation, this is a highly vectorized loop
        2539             :     !$acc parallel loop gang vector collapse(2) default(present)
        2540  9640925424 :     do k = 3, nz-1
        2541 >15915*10^7 :       do i = 1, ngrdcol
        2542 >29903*10^7 :         rhs(i,k) = 0.5_core_rknd  &
        2543 >14951*10^7 :                    * ( - lhs_diff(3,i,k) * xm(i,k-1)      &
        2544             :                        - lhs_diff(2,i,k) * xm(i,k)        &
        2545 >14951*10^7 :                        - lhs_diff(1,i,k) * xm(i,k+1) )    &
        2546             :                    + xm_tndcy(i,k)                        & ! RHS forcings
        2547 >75710*10^7 :                    + invrs_dt * xm(i,k)                     ! RHS time tendency
        2548             :       end do
        2549             :     end do
        2550             :     !$acc end parallel loop
        2551             :
        2552             :     ! Upper boundary calculation
        2553             :     !$acc parallel loop gang vector default(present)
        2554  1939530528 :     do i = 1, ngrdcol
        2555  3646749600 :       rhs(i,nz) = 0.5_core_rknd  &
        2556  1823374800 :                   * ( - lhs_diff(3,i,nz) * xm(i,nz-1)  &
        2557             :                       - lhs_diff(2,i,nz) * xm(i,nz) )  &
        2558             :                   + xm_tndcy(i,nz)                     & ! RHS forcings
        2559  5586280128 :                   + invrs_dt * xm(i,nz)                  ! RHS time tendency
        2560             :     end do
        2561             :     !$acc end parallel loop
        2562             :
        2563   116155728 :     if ( stats_metadata%l_stats_samp .and. ixm_ta > 0 ) then
        2564             :
        2565             :       !$acc update host( lhs_diff, xm )
        2566             :
        2567           0 :       do i = 1, ngrdcol
        2568             :
        2569             :         ! Statistics:  explicit contributions for um or vm.
        2570             :
        2571             :         ! xm term ta has both implicit and explicit components; call
        2572             :         ! stat_begin_update_pt.  Since stat_begin_update_pt automatically
        2573             :         ! subtracts the value sent in, reverse the sign on right-hand side
        2574             :         ! turbulent advection component.
        2575             :
        2576             :         ! Lower boundary
        2577             :         call stat_begin_update_pt( ixm_ta, 2, &                         ! intent(in)
        2578             :                                    0.5_core_rknd  &
        2579           0 :                                  * ( lhs_diff(2,i,2) * xm(i,2)   &      ! intent(in)
        2580           0 :                                  +   lhs_diff(1,i,2) * xm(i,3) ), &
        2581           0 :                                      stats_zt(i) )                      ! intent(inout)
        2582             :
        2583           0 :         do k = 3, nz-1
        2584             :
        2585             :           call stat_begin_update_pt( ixm_ta, k, &                         ! intent(in)
        2586             :                                      0.5_core_rknd  &
        2587           0 :                                    * ( lhs_diff(3,i,k) * xm(i,k-1) &
        2588             :                                    +   lhs_diff(2,i,k) * xm(i,k)   &      ! intent(in)
        2589           0 :                                    +   lhs_diff(1,i,k) * xm(i,k+1) ), &
        2590           0 :                                        stats_zt(i) )                      ! intent(inout)
        2591             :         end do
        2592             :
        2593             :         ! Upper boundary
        2594             :         call stat_begin_update_pt( ixm_ta, nz, &
        2595             :                                    0.5_core_rknd  &                       ! intent(in)
        2596           0 :                                  * ( lhs_diff(3,i,nz) * xm(i,nz-1) &
        2597             :                                  +   lhs_diff(2,i,nz) * xm(i,nz) ), &     ! intent(in)
        2598           0 :                                      stats_zt(i) )                        ! intent(inout)
        2599             :       end do
        2600             :     endif
        2601             :
        2602   116155728 :     if ( .not. l_imp_sfc_momentum_flux ) then
        2603             :
        2604             :       ! RHS generalized surface flux.
        2605             :       !$acc parallel loop gang vector default(present)
        2606  1939530528 :       do i = 1, ngrdcol
        2607  3646749600 :         rhs(i,2) = rhs(i,2) + invrs_rho_ds_zt(i,2)  &
        2608           0 :                               * gr%invrs_dzt(i,2)  &
        2609  3762905328 :                               * rho_ds_zm(i,1) * xpwp_sfc(i)
        2610             :       end do
        2611             :       !$acc end parallel loop
        2612             :
        2613   116155728 :       if ( stats_metadata%l_stats_samp .and. ixm_ta > 0 ) then
        2614             :
        2615             :         !$acc update host( invrs_rho_ds_zt, rho_ds_zm, xpwp_sfc )
        2616             :
        2617           0 :         do i = 1, ngrdcol
        2618             :           call stat_modify_pt( ixm_ta, 2,  &                    ! intent(in)
        2619           0 :                              + invrs_rho_ds_zt(i,2)  &
        2620           0 :                              * gr%invrs_dzt(i,2)  &
        2621           0 :                              * rho_ds_zm(i,1) * xpwp_sfc(i),  & ! intent(in)
        2622           0 :                                stats_zt(i) )                    ! intent(inout)
        2623             :         end do
        2624             :       end if
        2625             :
        2626             :     endif ! l_imp_sfc_momentum_flux
        2627             :
        2628   116155728 :     return
        2629             :
        2630             :   end subroutine windm_edsclrm_rhs
        2631             :
        2632             : end module advance_windm_edsclrm_module

Generated by: LCOV version 1.14