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