LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - clip_explicit.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 136 | 219 | 62.1 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 5 | 6 | 83.3 %

              Line data    Source code
           1             : !-------------------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module clip_explicit
           5             :
           6             :   implicit none
           7             :
           8             :   private
           9             :
          10             :   public :: clip_covars_denom, &
          11             :             clip_covar, &
          12             :             clip_covar_level, &
          13             :             clip_variance, &
          14             :             clip_skewness, &
          15             :             clip_skewness_core
          16             :
          17             :   ! Named constants to avoid string comparisons
          18             :   integer, parameter, public :: &
          19             :     clip_rtp2 = 1, &         ! Named constant for rtp2 clipping
          20             :     clip_thlp2 = 2, &        ! Named constant for thlp2 clipping
          21             :     clip_rtpthlp = 3, &      ! Named constant for rtpthlp clipping
          22             :     clip_up2 = 5, &          ! Named constant for up2 clipping
          23             :     clip_vp2 = 6, &          ! Named constant for vp2 clipping
          24             : !    clip_scalar = 7, &       ! Named constant for scalar clipping
          25             :     clip_wprtp = 8, &        ! Named constant for wprtp clipping
          26             :     clip_wpthlp = 9, &       ! Named constant for wpthlp clipping
          27             :     clip_upwp = 10, &        ! Named constant for upwp clipping
          28             :     clip_vpwp = 11, &        ! Named constant for vpwp clipping
          29             :     clip_wp2 = 12, &         ! Named constant for wp2 clipping
          30             :     clip_wpsclrp = 13, &     ! Named constant for wp scalar clipping
          31             :     clip_sclrp2 = 14, &      ! Named constant for sclrp2 clipping
          32             :     clip_sclrprtp = 15, &    ! Named constant for sclrprtp clipping
          33             :     clip_sclrpthlp = 16, &   ! Named constant for sclrpthlp clipping
          34             :     clip_wphydrometp = 17    ! Named constant for wphydrometp clipping
          35             :
          36             :   contains
          37             :
          38             :   !=============================================================================
          39    17870112 :   subroutine clip_covars_denom( nz, ngrdcol, sclr_dim, gr, dt, &
          40    17870112 :                                 rtp2, thlp2, up2, vp2, wp2, &
          41    17870112 :                                 sclrp2, wprtp_cl_num, wpthlp_cl_num, &
          42             :                                 wpsclrp_cl_num, upwp_cl_num, vpwp_cl_num, &
          43             :                                 l_predict_upwp_vpwp, &
          44             :                                 l_tke_aniso, &
          45             :                                 l_linearize_pbl_winds, &
          46             :                                 stats_metadata, &
          47    17870112 :                                 stats_zm, &
          48    17870112 :                                 wprtp, wpthlp, upwp, vpwp, wpsclrp, &
          49    17870112 :                                 upwp_pert, vpwp_pert )
          50             :
          51             :     ! Description:
          52             :     ! Some of the covariances found in the CLUBB model code need to be clipped
          53             :     ! multiple times during each timestep to ensure that the correlation between
          54             :     ! the two relevant variables stays between -1 and 1 at all times during the
          55             :     ! model run.  The covariances that need to be clipped multiple times are
          56             :     ! w'r_t', w'th_l', w'sclr', u'w', and v'w'.  One of the times that each one
          57             :     ! of these covariances is clipped is immediately after each one is set.
          58             :     ! However, each covariance still needs to be clipped two more times during
          59             :     ! each timestep (once after advance_xp2_xpyp is called and once after
          60             :     ! advance_wp2_wp3 is called).  This subroutine handles the times that the
          61             :     ! covariances are clipped away from the time that they are set.  In other
          62             :     ! words, this subroutine clips the covariances after the denominator terms
          63             :     ! in the relevant correlation equation have been altered, ensuring that
          64             :     ! all correlations will remain between -1 and 1 at all times.
          65             :
          66             :     ! References:
          67             :     ! None
          68             :     !-----------------------------------------------------------------------
          69             :
          70             :     use grid_class, only: &
          71             :         grid ! Type
          72             :
          73             :     use clubb_precision, only: &
          74             :         core_rknd ! Variable(s)
          75             :
          76             :     use stats_type, only: &
          77             :         stats ! Type
          78             :
          79             :      use stats_variables, only: &
          80             :         stats_metadata_type
          81             :
          82             :     implicit none
          83             :
          84             :     ! --------------------- Input Variables ---------------------
          85             :     integer, intent(in) :: &
          86             :       nz, &
          87             :       ngrdcol, &
          88             :       sclr_dim
          89             :
          90             :     type (grid), target, intent(in) :: gr
          91             :
          92             :     real( kind = core_rknd ), intent(in) :: &
          93             :       dt ! Timestep [s]
          94             :
          95             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
          96             :       rtp2,  & ! r_t'^2         [(kg/kg)^2]
          97             :       thlp2, & ! theta_l'^2     [K^2]
          98             :       up2,   & ! u'^2           [m^2/s^2]
          99             :       vp2,   & ! v'^2           [m^2/s^2]
         100             :       wp2      ! w'^2           [m^2/s^2]
         101             :
         102             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(in) :: &
         103             :       sclrp2 ! sclr'^2  [{units vary}^2]
         104             :
         105             :     integer, intent(in) :: &
         106             :       wprtp_cl_num,   &
         107             :       wpthlp_cl_num,  &
         108             :       wpsclrp_cl_num, &
         109             :       upwp_cl_num,    &
         110             :       vpwp_cl_num
         111             :
         112             :     logical, intent(in) :: &
         113             :       l_predict_upwp_vpwp,   & ! Flag to predict <u'w'> and <v'w'> along with <u> and <v> alongside
         114             :                                ! the advancement of <rt>, <w'rt'>, <thl>, <wpthlp>, <sclr>, and
         115             :                                ! <w'sclr'> in subroutine advance_xm_wpxp.  Otherwise, <u'w'> and
         116             :                                ! <v'w'> are still approximated by eddy diffusivity when <u> and <v>
         117             :                                ! are advanced in subroutine advance_windm_edsclrm.
         118             :       l_tke_aniso,           & ! For anisotropic turbulent kinetic energy, i.e. TKE = 1/2
         119             :                                ! (u'^2 + v'^2 + w'^2)
         120             :       l_linearize_pbl_winds    ! Flag (used by E3SM) to linearize PBL winds
         121             :
         122             :     type (stats_metadata_type), intent(in) :: &
         123             :       stats_metadata
         124             :
         125             :     ! --------------------- Input/Output Variables ---------------------
         126             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
         127             :       stats_zm
         128             :
         129             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
         130             :       wprtp,  & ! w'r_t'        [(kg/kg) m/s]
         131             :       wpthlp, & ! w'theta_l'    [K m/s]
         132             :       upwp,   & ! u'w'          [m^2/s^2]
         133             :       vpwp      ! v'w'          [m^2/s^2]
         134             :
         135             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(inout) :: &
         136             :       wpsclrp ! w'sclr'         [units m/s]
         137             :
         138             :     ! Variables used to track perturbed version of winds.
         139             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
         140             :       upwp_pert, & ! perturbed <u'w'> [m^2/s^2]
         141             :       vpwp_pert    ! perturbed <v'w'> [m^2/s^2]
         142             :
         143             :     ! --------------------- Local Variables ---------------------
         144             :     logical :: &
         145             :       l_first_clip_ts, & ! First instance of clipping in a timestep.
         146             :       l_last_clip_ts     ! Last instance of clipping in a timestep.
         147             :
         148             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         149    35740224 :       wprtp_chnge,  & ! Net change in w'r_t' due to clipping  [(kg/kg) m/s]
         150    35740224 :       wpthlp_chnge, & ! Net change in w'th_l' due to clipping [K m/s]
         151    35740224 :       upwp_chnge,   & ! Net change in u'w' due to clipping    [m^2/s^2]
         152    35740224 :       vpwp_chnge      ! Net change in v'w' due to clipping    [m^2/s^2]
         153             :
         154             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim) :: &
         155    35740224 :       wpsclrp_chnge   ! Net change in w'sclr' due to clipping [{units vary}]
         156             :
         157             :     integer :: sclr, i  ! scalar array index.
         158             :
         159             :     ! --------------------- Begin Code ---------------------
         160             :
         161             :     !$acc enter data create( wprtp_chnge, wpthlp_chnge, upwp_chnge, vpwp_chnge )
         162             :     !$acc enter data if( sclr_dim > 0 ) create( wpsclrp_chnge )
         163             :
         164             :     !!! Clipping for w'r_t'
         165             :     !
         166             :     ! Clipping w'r_t' at each vertical level, based on the
         167             :     ! correlation of w and r_t at each vertical level, such that:
         168             :     ! corr_(w,r_t) = w'r_t' / [ sqrt(w'^2) * sqrt(r_t'^2) ];
         169             :     ! -1 <= corr_(w,r_t) <= 1.
         170             :     !
         171             :     ! Since w'^2, r_t'^2, and w'r_t' are each advanced in different
         172             :     ! subroutines from each other in advance_clubb_core, clipping for w'r_t'
         173             :     ! is done three times during each timestep (once after each variable has
         174             :     ! been updated).
         175             :     !
         176             :     ! This subroutine handles the first and third instances of
         177             :     ! w'r_t' clipping.
         178             :     ! The first instance of w'r_t' clipping takes place after
         179             :     ! r_t'^2 is updated in advance_xp2_xpyp.
         180             :     ! The third instance of w'r_t' clipping takes place after
         181             :     ! w'^2 is updated in advance_wp2_wp3.
         182             :
         183             :     ! Used within subroutine clip_covar.
         184    17870112 :     if ( wprtp_cl_num == 1 ) then
         185           0 :       l_first_clip_ts = .true.
         186           0 :       l_last_clip_ts  = .false.
         187    17870112 :     elseif ( wprtp_cl_num == 2 ) then
         188     8935056 :       l_first_clip_ts = .false.
         189     8935056 :       l_last_clip_ts  = .false.
         190     8935056 :     elseif ( wprtp_cl_num == 3 ) then
         191     8935056 :       l_first_clip_ts = .false.
         192     8935056 :       l_last_clip_ts  = .true.
         193             :     endif
         194             :
         195             :     ! Clip w'r_t'
         196             :     call clip_covar( nz, ngrdcol, gr, clip_wprtp, l_first_clip_ts,  & ! intent(in)
         197             :                      l_last_clip_ts, dt, wp2, rtp2,                 & ! intent(in)
         198             :                      l_predict_upwp_vpwp,                           & ! intent(in)
         199             :                      stats_metadata,                                & ! intent(in)
         200             :                      stats_zm,                                      & ! intent(inout)
         201    17870112 :                      wprtp, wprtp_chnge )                             ! intent(inout)
         202             :
         203             :     !!! Clipping for w'th_l'
         204             :     !
         205             :     ! Clipping w'th_l' at each vertical level, based on the
         206             :     ! correlation of w and th_l at each vertical level, such that:
         207             :     ! corr_(w,th_l) = w'th_l' / [ sqrt(w'^2) * sqrt(th_l'^2) ];
         208             :     ! -1 <= corr_(w,th_l) <= 1.
         209             :     !
         210             :     ! Since w'^2, th_l'^2, and w'th_l' are each advanced in different
         211             :     ! subroutines from each other in advance_clubb_core, clipping for w'th_l'
         212             :     ! is done three times during each timestep (once after each variable has
         213             :     ! been updated).
         214             :     !
         215             :     ! This subroutine handles the first and third instances of
         216             :     ! w'th_l' clipping.
         217             :     ! The first instance of w'th_l' clipping takes place after
         218             :     ! th_l'^2 is updated in advance_xp2_xpyp.
         219             :     ! The third instance of w'th_l' clipping takes place after
         220             :     ! w'^2 is updated in advance_wp2_wp3.
         221             :
         222             :     ! Used within subroutine clip_covar.
         223    17870112 :     if ( wpthlp_cl_num == 1 ) then
         224           0 :       l_first_clip_ts = .true.
         225           0 :       l_last_clip_ts  = .false.
         226    17870112 :     elseif ( wpthlp_cl_num == 2 ) then
         227     8935056 :       l_first_clip_ts = .false.
         228     8935056 :       l_last_clip_ts  = .false.
         229     8935056 :     elseif ( wpthlp_cl_num == 3 ) then
         230     8935056 :       l_first_clip_ts = .false.
         231     8935056 :       l_last_clip_ts  = .true.
         232             :     endif
         233             :
         234             :     ! Clip w'th_l'
         235             :     call clip_covar( nz, ngrdcol, gr, clip_wpthlp, l_first_clip_ts, & ! intent(in)
         236             :                      l_last_clip_ts, dt, wp2, thlp2,                & ! intent(in)
         237             :                      l_predict_upwp_vpwp,                           & ! intent(in)
         238             :                      stats_metadata,                                & ! intent(in)
         239             :                      stats_zm,                                      & ! intent(inout)
         240    17870112 :                      wpthlp, wpthlp_chnge )                           ! intent(inout)
         241             :
         242             :     !!! Clipping for w'sclr'
         243             :     !
         244             :     ! Clipping w'sclr' at each vertical level, based on the
         245             :     ! correlation of w and sclr at each vertical level, such that:
         246             :     ! corr_(w,sclr) = w'sclr' / [ sqrt(w'^2) * sqrt(sclr'^2) ];
         247             :     ! -1 <= corr_(w,sclr) <= 1.
         248             :     !
         249             :     ! Since w'^2, sclr'^2, and w'sclr' are each advanced in different
         250             :     ! subroutines from each other in advance_clubb_core, clipping for w'sclr'
         251             :     ! is done three times during each timestep (once after each variable has
         252             :     ! been updated).
         253             :     !
         254             :     ! This subroutine handles the first and third instances of
         255             :     ! w'sclr' clipping.
         256             :     ! The first instance of w'sclr' clipping takes place after
         257             :     ! sclr'^2 is updated in advance_xp2_xpyp.
         258             :     ! The third instance of w'sclr' clipping takes place after
         259             :     ! w'^2 is updated in advance_wp2_wp3.
         260             :
         261             :     ! Used within subroutine clip_covar.
         262    17870112 :     if ( wpsclrp_cl_num == 1 ) then
         263           0 :       l_first_clip_ts = .true.
         264           0 :       l_last_clip_ts  = .false.
         265    17870112 :     elseif ( wpsclrp_cl_num == 2 ) then
         266     8935056 :       l_first_clip_ts = .false.
         267     8935056 :       l_last_clip_ts  = .false.
         268     8935056 :     elseif ( wpsclrp_cl_num == 3 ) then
         269     8935056 :       l_first_clip_ts = .false.
         270     8935056 :       l_last_clip_ts  = .true.
         271             :     endif
         272             :
         273             :     ! Clip w'sclr'
         274    17870112 :     do sclr = 1, sclr_dim
         275             :       call clip_covar( nz, ngrdcol, gr, clip_wpsclrp, l_first_clip_ts,  & ! intent(in)
         276             :                        l_last_clip_ts, dt, wp2(:,:), sclrp2(:,:,sclr),  & ! intent(in)
         277             :                        l_predict_upwp_vpwp,                             & ! intent(in)
         278             :                        stats_metadata,                                  & ! intent(in)
         279             :                        stats_zm,                                        & ! intent(inout)
         280    17870112 :                        wpsclrp(:,:,sclr), wpsclrp_chnge(:,:,sclr) )       ! intent(inout)
         281             :     enddo
         282             :
         283             :
         284             :     !!! Clipping for u'w'
         285             :     !
         286             :     ! Clipping u'w' at each vertical level, based on the
         287             :     ! correlation of u and w at each vertical level, such that:
         288             :     ! corr_(u,w) = u'w' / [ sqrt(u'^2) * sqrt(w'^2) ];
         289             :     ! -1 <= corr_(u,w) <= 1.
         290             :     !
         291             :     ! Since w'^2, u'^2, and u'w' are each advanced in different
         292             :     ! subroutines from each other in advance_clubb_core, clipping for u'w'
         293             :     ! is done three times during each timestep (once after each variable has
         294             :     ! been updated).
         295             :     !
         296             :     ! This subroutine handles the first and second instances of
         297             :     ! u'w' clipping.
         298             :     ! The first instance of u'w' clipping takes place after
         299             :     ! u'^2 is updated in advance_xp2_xpyp.
         300             :     ! The second instance of u'w' clipping takes place after
         301             :     ! w'^2 is updated in advance_wp2_wp3.
         302             :
         303             :     ! Used within subroutine clip_covar.
         304    17870112 :     if ( upwp_cl_num == 1 ) then
         305           0 :       l_first_clip_ts = .true.
         306           0 :       l_last_clip_ts  = .false.
         307    17870112 :     elseif ( upwp_cl_num == 2 ) then
         308     8935056 :       l_first_clip_ts = .false.
         309     8935056 :       l_last_clip_ts  = .false.
         310     8935056 :     elseif ( upwp_cl_num == 3 ) then
         311     8935056 :       l_first_clip_ts = .false.
         312     8935056 :       l_last_clip_ts  = .true.
         313             :     endif
         314             :
         315             :     ! Clip u'w'
         316    17870112 :     if ( l_tke_aniso ) then
         317             :       call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
         318             :                        l_last_clip_ts, dt, wp2, up2,                & ! intent(in)
         319             :                        l_predict_upwp_vpwp,                         & ! intent(in)
         320             :                        stats_metadata,                              & ! intent(in)
         321             :                        stats_zm,                                    & ! intent(inout)
         322    17870112 :                        upwp, upwp_chnge )                             ! intent(inout)
         323             :
         324    17870112 :       if ( l_linearize_pbl_winds ) then
         325             :         call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
         326             :                          l_last_clip_ts, dt, wp2, up2,                & ! intent(in)
         327             :                          l_predict_upwp_vpwp,                         & ! intent(in)
         328             :                          stats_metadata,                              & ! intent(in)
         329             :                          stats_zm,                                    & ! intent(inout)
         330           0 :                          upwp_pert, upwp_chnge )                        ! intent(inout)
         331             :       endif ! l_linearize_pbl_winds
         332             :     else
         333             :       ! In this case, up2 = wp2, and the variable `up2' does not interact
         334             :       call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
         335             :                        l_last_clip_ts, dt, wp2, wp2,                & ! intent(in)
         336             :                        l_predict_upwp_vpwp,                         & ! intent(in)
         337             :                        stats_metadata,                              & ! intent(in)
         338             :                        stats_zm,                                    & ! intent(inout)
         339           0 :                        upwp, upwp_chnge )                             ! intent(inout)
         340             :
         341           0 :       if ( l_linearize_pbl_winds ) then
         342             :           call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
         343             :                            l_last_clip_ts, dt, wp2, wp2,                & ! intent(in)
         344             :                            l_predict_upwp_vpwp,                         & ! intent(in)
         345             :                            stats_metadata,                              & ! intent(in)
         346             :                            stats_zm,                                    & ! intent(inout)
         347           0 :                            upwp_pert, upwp_chnge )                        ! intent(inout)
         348             :       endif ! l_linearize_pbl_winds
         349             :     end if
         350             :
         351             :
         352             :
         353             :     !!! Clipping for v'w'
         354             :     !
         355             :     ! Clipping v'w' at each vertical level, based on the
         356             :     ! correlation of v and w at each vertical level, such that:
         357             :     ! corr_(v,w) = v'w' / [ sqrt(v'^2) * sqrt(w'^2) ];
         358             :     ! -1 <= corr_(v,w) <= 1.
         359             :     !
         360             :     ! Since w'^2, v'^2, and v'w' are each advanced in different
         361             :     ! subroutines from each other in advance_clubb_core, clipping for v'w'
         362             :     ! is done three times during each timestep (once after each variable has
         363             :     ! been updated).
         364             :     !
         365             :     ! This subroutine handles the first and second instances of
         366             :     ! v'w' clipping.
         367             :     ! The first instance of v'w' clipping takes place after
         368             :     ! v'^2 is updated in advance_xp2_xpyp.
         369             :     ! The second instance of v'w' clipping takes place after
         370             :     ! w'^2 is updated in advance_wp2_wp3.
         371             :
         372             :     ! Used within subroutine clip_covar.
         373    17870112 :     if ( vpwp_cl_num == 1 ) then
         374           0 :       l_first_clip_ts = .true.
         375           0 :       l_last_clip_ts  = .false.
         376    17870112 :     elseif ( vpwp_cl_num == 2 ) then
         377     8935056 :       l_first_clip_ts = .false.
         378     8935056 :       l_last_clip_ts  = .false.
         379     8935056 :     elseif ( vpwp_cl_num == 3 ) then
         380     8935056 :       l_first_clip_ts = .false.
         381     8935056 :       l_last_clip_ts  = .true.
         382             :     endif
         383             :
         384    17870112 :     if ( l_tke_aniso ) then
         385             :       call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
         386             :                        l_last_clip_ts, dt, wp2, vp2,                & ! intent(in)
         387             :                        l_predict_upwp_vpwp,                         & ! intent(in)
         388             :                        stats_metadata,                              & ! intent(in)
         389             :                        stats_zm,                                    & ! intent(inout)
         390    17870112 :                        vpwp, vpwp_chnge )                             ! intent(inout)
         391             :
         392    17870112 :       if ( l_linearize_pbl_winds ) then
         393             :         call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
         394             :                          l_last_clip_ts, dt, wp2, vp2,                & ! intent(in)
         395             :                          l_predict_upwp_vpwp,                         & ! intent(in)
         396             :                          stats_metadata,                              & ! intent(in)
         397             :                          stats_zm,                                    & ! intent(inout)
         398           0 :                          vpwp_pert, vpwp_chnge )                        ! intent(inout)
         399             :       endif ! l_linearize_pbl_winds
         400             :     else
         401             :       ! In this case, vp2 = wp2, and the variable `vp2' does not interact
         402             :       call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
         403             :                        l_last_clip_ts, dt, wp2, wp2,                & ! intent(in)
         404             :                        l_predict_upwp_vpwp,                         & ! intent(in)
         405             :                        stats_metadata,                              & ! intent(in)
         406             :                        stats_zm,                                    & ! intent(inout)
         407           0 :                        vpwp, vpwp_chnge )                             ! intent(inout)
         408             :
         409           0 :       if ( l_linearize_pbl_winds ) then
         410             :         call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
         411             :                          l_last_clip_ts, dt, wp2, wp2,                & ! intent(in)
         412             :                          l_predict_upwp_vpwp,                         & ! intent(in)
         413             :                          stats_metadata,                              & ! intent(in)                       stats_metadata,                              & ! intent(in)
         414             :                          stats_zm,                                    & ! intent(inout)
         415           0 :                          vpwp_pert, vpwp_chnge )                        ! intent(inout)
         416             :       endif ! l_linearize_pbl_winds
         417             :     end if
         418             :
         419             :     !$acc exit data delete( wprtp_chnge, wpthlp_chnge, upwp_chnge, vpwp_chnge )
         420             :     !$acc exit data if( sclr_dim > 0 ) delete( wpsclrp_chnge )
         421             :
         422    17870112 :     return
         423             :   end subroutine clip_covars_denom
         424             :
         425             :   !=============================================================================
         426   116155728 :   subroutine clip_covar( nz, ngrdcol, gr, solve_type, l_first_clip_ts,  &
         427   116155728 :                          l_last_clip_ts, dt, xp2, yp2,  &
         428             :                          l_predict_upwp_vpwp, &
         429             :                          stats_metadata, &
         430   116155728 :                          stats_zm, &
         431   116155728 :                          xpyp, xpyp_chnge )
         432             :
         433             :     ! Description:
         434             :     ! Clipping the value of covariance x'y' based on the correlation between x
         435             :     ! and y.
         436             :     !
         437             :     ! The correlation between variables x and y is:
         438             :     !
         439             :     ! corr_(x,y) = x'y' / [ sqrt(x'^2) * sqrt(y'^2) ];
         440             :     !
         441             :     ! where x'^2 is the variance of x, y'^2 is the variance of y, and x'y' is
         442             :     ! the covariance of x and y.
         443             :     !
         444             :     ! The correlation of two variables must always have a value between -1
         445             :     ! and 1, such that:
         446             :     !
         447             :     ! -1 <= corr_(x,y) <= 1.
         448             :     !
         449             :     ! Therefore, there is an upper limit on x'y', such that:
         450             :     !
         451             :     ! x'y' <=  [ sqrt(x'^2) * sqrt(y'^2) ];
         452             :     !
         453             :     ! and a lower limit on x'y', such that:
         454             :     !
         455             :     ! x'y' >= -[ sqrt(x'^2) * sqrt(y'^2) ].
         456             :     !
         457             :     ! The values of x'y', x'^2, and y'^2 are all found on momentum levels.
         458             :     !
         459             :     ! The value of x'y' may need to be clipped whenever x'y', x'^2, or y'^2 is
         460             :     ! updated.
         461             :     !
         462             :     ! The following covariances are found in the code:
         463             :     !
         464             :     ! w'r_t', w'th_l', w'sclr', (computed in advance_xm_wpxp);
         465             :     ! r_t'th_l', sclr'r_t', sclr'th_l', (computed in advance_xp2_xpyp);
         466             :     ! u'w', v'w', w'edsclr' (computed in advance_windm_edsclrm);
         467             :     ! and w'hm' (computed in setup_pdf_parameters).
         468             :
         469             :     ! References:
         470             :     ! None
         471             :     !-----------------------------------------------------------------------
         472             :
         473             :     use grid_class, only: &
         474             :         grid ! Type
         475             :
         476             :     use constants_clubb, only: &
         477             :         max_mag_correlation,      & ! Constant(s)
         478             :         max_mag_correlation_flux
         479             :
         480             :     use clubb_precision, only: &
         481             :         core_rknd ! Variable(s)
         482             :
         483             :     use stats_type_utilities, only: &
         484             :         stat_begin_update,  & ! Procedure(s)
         485             :         stat_modify, &
         486             :         stat_end_update
         487             :
         488             :     use stats_variables, only: &
         489             :         stats_metadata_type
         490             :
         491             :     use stats_type, only: stats ! Type
         492             :
         493             :     implicit none
         494             :
         495             :     ! -------------------------- Input Variables --------------------------
         496             :     integer, intent(in) :: &
         497             :       nz, &
         498             :       ngrdcol
         499             :
         500             :     type (grid), target, intent(in) :: gr
         501             :
         502             :     integer, intent(in) :: &
         503             :       solve_type       ! Variable being solved; used for STATS.
         504             :
         505             :     logical, intent(in) :: &
         506             :       l_first_clip_ts, & ! First instance of clipping in a timestep.
         507             :       l_last_clip_ts     ! Last instance of clipping in a timestep.
         508             :
         509             :     real( kind = core_rknd ), intent(in) ::  &
         510             :       dt     ! Model timestep; used here for STATS           [s]
         511             :
         512             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         513             :       xp2, & ! Variance of x, x'^2 (momentum levels)         [{x units}^2]
         514             :       yp2    ! Variance of y, y'^2 (momentum levels)         [{y units}^2]
         515             :
         516             :     logical, intent(in) :: &
         517             :       l_predict_upwp_vpwp ! Flag to predict <u'w'> and <v'w'> along with <u> and <v> alongside the
         518             :                           ! advancement of <rt>, <w'rt'>, <thl>, <wpthlp>, <sclr>, and <w'sclr'> in
         519             :                           ! subroutine advance_xm_wpxp.  Otherwise, <u'w'> and <v'w'> are still
         520             :                           ! approximated by eddy diffusivity when <u> and <v> are advanced in
         521             :                           ! subroutine advance_windm_edsclrm.
         522             :
         523             :     type (stats_metadata_type), intent(in) :: &
         524             :       stats_metadata
         525             :
         526             :     ! -------------------------- InOut Variables --------------------------
         527             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
         528             :       stats_zm
         529             :
         530             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
         531             :       xpyp   ! Covariance of x and y, x'y' (momentum levels) [{x units}*{y units}]
         532             :
         533             :     !-------------------------- Output Variable --------------------------
         534             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         535             :       xpyp_chnge  ! Net change in x'y' due to clipping [{x units}*{y units}]
         536             :
         537             :
         538             :     ! -------------------------- Local Variables --------------------------
         539             :     real( kind = core_rknd ) ::  &
         540             :       max_mag_corr, &    ! Maximum magnitude of a correlation allowed
         541             :       xpyp_bound
         542             :
         543             :     integer :: i, k  ! Array index
         544             :
         545             :     integer :: &
         546             :       ixpyp_cl
         547             :
         548             :     ! -------------------------- Begin Code --------------------------
         549             :
         550   142960896 :     select case ( solve_type )
         551             :     case ( clip_wprtp )   ! wprtp clipping budget term
         552    26805168 :       ixpyp_cl = stats_metadata%iwprtp_cl
         553             :     case ( clip_wpthlp )   ! wpthlp clipping budget term
         554    26805168 :       ixpyp_cl = stats_metadata%iwpthlp_cl
         555             :     case ( clip_rtpthlp )   ! rtpthlp clipping budget term
         556     8935056 :       ixpyp_cl = stats_metadata%irtpthlp_cl
         557             :     case ( clip_upwp )   ! upwp clipping budget term
         558    26805168 :       if ( l_predict_upwp_vpwp ) then
         559    26805168 :         ixpyp_cl = stats_metadata%iupwp_cl
         560             :       else
         561           0 :         ixpyp_cl = 0
         562             :       endif ! l_predict_upwp_vpwp
         563             :     case ( clip_vpwp )   ! vpwp clipping budget term
         564    26805168 :       if ( l_predict_upwp_vpwp ) then
         565    26805168 :         ixpyp_cl = stats_metadata%ivpwp_cl
         566             :       else
         567           0 :         ixpyp_cl = 0
         568             :       endif ! l_predict_upwp_vpwp
         569             :     case default   ! scalars (or upwp/vpwp) are involved
         570   116155728 :       ixpyp_cl = 0
         571             :     end select
         572             :
         573             :
         574   116155728 :     if ( stats_metadata%l_stats_samp ) then
         575             :
         576             :       !$acc update host( xpyp )
         577             :
         578           0 :       if ( l_first_clip_ts ) then
         579           0 :         do i = 1, ngrdcol
         580           0 :           call stat_begin_update( nz, ixpyp_cl, xpyp(i,:) / dt, & ! intent(in)
         581           0 :                                   stats_zm(i) ) ! intent(inout)
         582             :         end do
         583             :       else
         584           0 :         do i = 1, ngrdcol
         585           0 :           call stat_modify( nz, ixpyp_cl, -xpyp(i,:) / dt, & ! intent(in)
         586           0 :                             stats_zm(i) ) ! intent(inout)
         587             :         end do
         588             :       endif
         589             :     endif
         590             :
         591             :     ! When clipping for wprtp or wpthlp, use the special value for
         592             :     ! max_mag_correlation_flux.  For all other correlations, use
         593             :     ! max_mag_correlation.
         594             :     if ( ( solve_type == clip_wprtp ) .or. ( solve_type == clip_wpthlp ) ) then
         595             :        max_mag_corr = max_mag_correlation_flux
         596             :     else ! All other covariances
         597             :        max_mag_corr = max_mag_correlation
         598             :     endif ! solve_type
         599             :
         600             :     ! The value of x'y' at the surface (or lower boundary) is a set value that
         601             :     ! is either specified or determined elsewhere in a surface subroutine.  It
         602             :     ! is ensured elsewhere that the correlation between x and y at the surface
         603             :     ! (or lower boundary) is between -1 and 1.  Thus, the covariance clipping
         604             :     ! code does not need to be invoked at the lower boundary.  Likewise, the
         605             :     ! value of x'y' is set at the upper boundary, so the covariance clipping
         606             :     ! code does not need to be invoked at the upper boundary.
         607             :     ! Note that if clipping were applied at the lower boundary, momentum will
         608             :     ! not be conserved, therefore it should never be added.
         609             :     !$acc parallel loop gang vector collapse(2) default(present)
         610  9757081152 :     do k = 2, nz-1
         611 >16109*10^7 :       do i = 1, ngrdcol
         612 >15134*10^7 :         xpyp_bound = max_mag_corr * sqrt( xp2(i,k) * yp2(i,k) )
         613             :
         614             :         ! Clipping for xpyp at an upper limit corresponding with a correlation
         615             :         ! between x and y of max_mag_corr.
         616 >16098*10^7 :         if ( xpyp(i,k) > xpyp_bound ) then
         617             :
         618    38138270 :           xpyp_chnge(i,k) = xpyp_bound - xpyp(i,k)
         619    38138270 :           xpyp(i,k) = xpyp_bound
         620             :
         621             :         ! Clipping for xpyp at a lower limit corresponding with a correlation
         622             :         ! between x and y of -max_mag_corr.
         623 >15130*10^7 :         else if ( xpyp(i,k) < -xpyp_bound ) then
         624             :
         625    34845712 :           xpyp_chnge(i,k) = -xpyp_bound - xpyp(i,k)
         626    34845712 :           xpyp(i,k) = -xpyp_bound
         627             :
         628             :         else
         629             :
         630 >15126*10^7 :           xpyp_chnge(i,k) = 0.0_core_rknd
         631             :
         632             :         end if
         633             :       end do
         634             :     end do
         635             :     !$acc end parallel loop
         636             :
         637             :     ! Since there is no covariance clipping at the upper or lower boundaries,
         638             :     ! the change in x'y' due to covariance clipping at those levels is 0.
         639             :     !$acc parallel loop gang vector default(present)
         640  1939530528 :     do i = 1, ngrdcol
         641  1823374800 :       xpyp_chnge(i,1)  = 0.0_core_rknd
         642  1939530528 :       xpyp_chnge(i,nz) = 0.0_core_rknd
         643             :     end do
         644             :     !$acc end parallel loop
         645             :
         646   116155728 :     if ( stats_metadata%l_stats_samp ) then
         647             :
         648             :       !$acc update host( xpyp )
         649             :
         650           0 :       if ( l_last_clip_ts ) then
         651           0 :         do i = 1, ngrdcol
         652           0 :           call stat_end_update( nz, ixpyp_cl, xpyp(i,:) / dt, & ! intent(in)
         653           0 :                                 stats_zm(i) ) ! intent(inout)
         654             :         end do
         655             :       else
         656           0 :         do i = 1, ngrdcol
         657           0 :           call stat_modify( nz, ixpyp_cl, xpyp(i,:) / dt, & ! intent(in)
         658           0 :                             stats_zm(i) ) ! intent(inout)
         659             :         end do
         660             :       endif
         661             :     endif
         662             :
         663   116155728 :     return
         664             :
         665             :   end subroutine clip_covar
         666             :
         667             :   !=============================================================================
         668           0 :   subroutine clip_covar_level( solve_type, level, l_first_clip_ts,  &
         669             :                                l_last_clip_ts, dt, xp2, yp2,  &
         670             :                                l_predict_upwp_vpwp, &
         671             :                                stats_metadata, &
         672             :                                stats_zm, &
         673             :                                xpyp, xpyp_chnge )
         674             :
         675             :     ! Description:
         676             :     ! Clipping the value of covariance x'y' based on the correlation between x
         677             :     ! and y.  This is all done at a single vertical level.
         678             :     !
         679             :     ! The correlation between variables x and y is:
         680             :     !
         681             :     ! corr_(x,y) = x'y' / [ sqrt(x'^2) * sqrt(y'^2) ];
         682             :     !
         683             :     ! where x'^2 is the variance of x, y'^2 is the variance of y, and x'y' is
         684             :     ! the covariance of x and y.
         685             :     !
         686             :     ! The correlation of two variables must always have a value between -1
         687             :     ! and 1, such that:
         688             :     !
         689             :     ! -1 <= corr_(x,y) <= 1.
         690             :     !
         691             :     ! Therefore, there is an upper limit on x'y', such that:
         692             :     !
         693             :     ! x'y' <=  [ sqrt(x'^2) * sqrt(y'^2) ];
         694             :     !
         695             :     ! and a lower limit on x'y', such that:
         696             :     !
         697             :     ! x'y' >= -[ sqrt(x'^2) * sqrt(y'^2) ].
         698             :     !
         699             :     ! The values of x'y', x'^2, and y'^2 are all found on momentum levels.
         700             :     !
         701             :     ! The value of x'y' may need to be clipped whenever x'y', x'^2, or y'^2 is
         702             :     ! updated.
         703             :     !
         704             :     ! The following covariances are found in the code:
         705             :     !
         706             :     ! w'r_t', w'th_l', w'sclr', (computed in advance_xm_wpxp);
         707             :     ! r_t'th_l', sclr'r_t', sclr'th_l', (computed in advance_xp2_xpyp);
         708             :     ! u'w', v'w', w'edsclr' (computed in advance_windm_edsclrm);
         709             :     ! and w'hm' (computed in setup_pdf_parameters).
         710             :
         711             :     ! References:
         712             :     ! None
         713             :     !-----------------------------------------------------------------------
         714             :
         715             :     use constants_clubb, only: &
         716             :         max_mag_correlation,      & ! Constant(s)
         717             :         max_mag_correlation_flux, &
         718             :         zero
         719             :
         720             :     use clubb_precision, only: &
         721             :         core_rknd ! Variable(s)
         722             :
         723             :     use stats_type_utilities, only: &
         724             :         stat_begin_update_pt, & ! Procedure(s)
         725             :         stat_modify_pt,       &
         726             :         stat_end_update_pt
         727             :
         728             :     use stats_variables, only: &
         729             :         stats_metadata_type
         730             :
         731             :     use stats_type, only: stats ! Type
         732             :
         733             :     implicit none
         734             :
         735             :     type (stats), target, intent(inout) :: &
         736             :       stats_zm
         737             :
         738             :     !------------------------- Input Variables -------------------------
         739             :     integer, intent(in) :: &
         740             :       solve_type, & ! Variable being solved; used for STATS
         741             :       level         ! Vertical level index
         742             :
         743             :     logical, intent(in) :: &
         744             :       l_first_clip_ts, & ! First instance of clipping in a timestep.
         745             :       l_last_clip_ts     ! Last instance of clipping in a timestep.
         746             :
         747             :     real( kind = core_rknd ), intent(in) ::  &
         748             :       dt     ! Model timestep; used here for STATS        [s]
         749             :
         750             :     real( kind = core_rknd ), intent(in) :: &
         751             :       xp2, & ! Variance of x, <x'^2>                      [{x units}^2]
         752             :       yp2    ! Variance of y, <y'^2>                      [{y units}^2]
         753             :
         754             :     logical, intent(in) :: &
         755             :       l_predict_upwp_vpwp ! Flag to predict <u'w'> and <v'w'> along with <u> and <v> alongside the
         756             :                           ! advancement of <rt>, <w'rt'>, <thl>, <wpthlp>, <sclr>, and <w'sclr'> in
         757             :                           ! subroutine advance_xm_wpxp.  Otherwise, <u'w'> and <v'w'> are still
         758             :                           ! approximated by eddy diffusivity when <u> and <v> are advanced in
         759             :                           ! subroutine advance_windm_edsclrm.
         760             :
         761             :     type (stats_metadata_type), intent(in) :: &
         762             :       stats_metadata
         763             :
         764             :     !------------------------- InOut Variable -------------------------
         765             :     real( kind = core_rknd ), intent(inout) :: &
         766             :       xpyp   ! Covariance of x and y, <x'y'>              [{x units}*{y units}]
         767             :
         768             :     !------------------------- Output Variable -------------------------
         769             :     real( kind = core_rknd ), intent(out) :: &
         770             :       xpyp_chnge  ! Net change in <x'y'> due to clipping  [{x units}*{y units}]
         771             :
         772             :
         773             :     !------------------------- Local Variables -------------------------
         774             :     real( kind = core_rknd ) ::  &
         775             :       max_mag_corr    ! Maximum magnitude of a correlation allowed
         776             :
         777             :     integer :: &
         778             :       ixpyp_cl    ! Statistics index
         779             :
         780             :     !------------------------- Begin Code -------------------------
         781             :
         782           0 :     select case ( solve_type )
         783             :     case ( clip_wprtp )   ! wprtp clipping budget term
         784           0 :       ixpyp_cl = stats_metadata%iwprtp_cl
         785             :     case ( clip_wpthlp )   ! wpthlp clipping budget term
         786           0 :       ixpyp_cl = stats_metadata%iwpthlp_cl
         787             :     case ( clip_rtpthlp )   ! rtpthlp clipping budget term
         788           0 :       ixpyp_cl = stats_metadata%irtpthlp_cl
         789             :     case ( clip_upwp )   ! upwp clipping budget term
         790           0 :       if ( l_predict_upwp_vpwp ) then
         791           0 :         ixpyp_cl = stats_metadata%iupwp_cl
         792             :       else
         793           0 :         ixpyp_cl = 0
         794             :       endif ! l_predict_upwp_vpwp
         795             :     case ( clip_vpwp )   ! vpwp clipping budget term
         796           0 :       if ( l_predict_upwp_vpwp ) then
         797           0 :         ixpyp_cl = stats_metadata%ivpwp_cl
         798             :       else
         799           0 :         ixpyp_cl = 0
         800             :       endif ! l_predict_upwp_vpwp
         801             :     case default   ! scalars (or upwp/vpwp) are involved
         802           0 :       ixpyp_cl = 0
         803             :     end select
         804             :
         805             :
         806           0 :     if ( stats_metadata%l_stats_samp ) then
         807           0 :        if ( l_first_clip_ts ) then
         808             :           call stat_begin_update_pt( ixpyp_cl, level, xpyp / dt, & ! intent(in)
         809           0 :                                      stats_zm ) ! intent(inout)
         810             :        else
         811             :           call stat_modify_pt( ixpyp_cl, level, -xpyp / dt, & ! intent(in)
         812           0 :                                stats_zm ) ! intent(inout)
         813             :        endif
         814             :     endif
         815             :
         816             :     ! When clipping for wprtp or wpthlp, use the special value for
         817             :     ! max_mag_correlation_flux.  For all other correlations, use
         818             :     ! max_mag_correlation.
         819             :     if ( ( solve_type == clip_wprtp ) .or. ( solve_type == clip_wpthlp ) ) then
         820             :        max_mag_corr = max_mag_correlation_flux
         821             :     else ! All other covariances
         822             :        max_mag_corr = max_mag_correlation
         823             :     endif ! solve_type
         824             :
         825             :     ! The value of x'y' at the surface (or lower boundary) is a set value that
         826             :     ! is either specified or determined elsewhere in a surface subroutine.  It
         827             :     ! is ensured elsewhere that the correlation between x and y at the surface
         828             :     ! (or lower boundary) is between -1 and 1.  Thus, the covariance clipping
         829             :     ! code does not need to be invoked at the lower boundary.  Likewise, the
         830             :     ! value of x'y' is set at the upper boundary, so the covariance clipping
         831             :     ! code does not need to be invoked at the upper boundary.
         832             :     ! Note that if clipping were applied at the lower boundary, momentum will
         833             :     ! not be conserved, therefore it should never be added.
         834             :
         835             :     ! Clipping for xpyp at an upper limit corresponding with a correlation
         836             :     ! between x and y of max_mag_corr.
         837           0 :     if ( xpyp > max_mag_corr * sqrt( xp2 * yp2 ) ) then
         838             :
         839           0 :         xpyp_chnge = max_mag_corr * sqrt( xp2 * yp2 ) - xpyp
         840             :
         841           0 :         xpyp = max_mag_corr * sqrt( xp2 * yp2 )
         842             :
         843             :     ! Clipping for xpyp at a lower limit corresponding with a correlation
         844             :     ! between x and y of -max_mag_corr.
         845           0 :     elseif ( xpyp < -max_mag_corr * sqrt( xp2 * yp2 ) ) then
         846             :
         847           0 :         xpyp_chnge = -max_mag_corr * sqrt( xp2 * yp2 ) - xpyp
         848             :
         849           0 :         xpyp = -max_mag_corr * sqrt( xp2 * yp2 )
         850             :
         851             :     else
         852             :
         853           0 :         xpyp_chnge = zero
         854             :
         855             :     endif
         856             :
         857           0 :     if ( stats_metadata%l_stats_samp ) then
         858           0 :        if ( l_last_clip_ts ) then
         859             :           call stat_end_update_pt( ixpyp_cl, level, xpyp / dt, & ! intent(in)
         860           0 :                                    stats_zm ) ! intent(inout)
         861             :        else
         862             :           call stat_modify_pt( ixpyp_cl, level, xpyp / dt, & ! intent(in)
         863           0 :                                stats_zm ) ! intent(inout)
         864             :        endif
         865             :     endif
         866             :
         867             :
         868           0 :     return
         869             :   end subroutine clip_covar_level
         870             :
         871             :   !=============================================================================
         872    44675280 :   subroutine clip_variance( nz, ngrdcol, gr, solve_type, dt, threshold, &
         873             :                             stats_metadata, &
         874    44675280 :                             stats_zm, &
         875    44675280 :                             xp2 )
         876             :
         877             :     ! Description:
         878             :     ! Clipping the value of variance x'^2 based on a minimum threshold value.
         879             :     ! The threshold value must be greater than or equal to 0.
         880             :     !
         881             :     ! The values of x'^2 are found on the momentum levels.
         882             :     !
         883             :     ! The following variances are found in the code:
         884             :     !
         885             :     ! r_t'^2, th_l'^2, u'^2, v'^2, sclr'^2, (computed in advance_xp2_xpyp);
         886             :     ! w'^2 (computed in advance_wp2_wp3).
         887             :
         888             :     ! References:
         889             :     ! None
         890             :     !-----------------------------------------------------------------------
         891             :
         892             :     use grid_class, only: &
         893             :         grid ! Type
         894             :
         895             :     use clubb_precision, only: &
         896             :         core_rknd ! Variable(s)
         897             :
         898             :     use stats_type_utilities, only: &
         899             :         stat_begin_update,  & ! Procedure(s)
         900             :         stat_end_update
         901             :
         902             :     use stats_variables, only: &
         903             :         stats_metadata_type
         904             :
         905             :     use stats_type, only: stats ! Type
         906             :
         907             :     implicit none
         908             :
         909             :     ! -------------------- Input Variables --------------------
         910             :     integer, intent(in) :: &
         911             :       nz, &
         912             :       ngrdcol
         913             :
         914             :     type (grid), target, intent(in) :: gr
         915             :
         916             :     integer, intent(in) :: &
         917             :       solve_type  ! Variable being solved; used for STATS.
         918             :
         919             :     real( kind = core_rknd ), intent(in) :: &
         920             :       dt          ! Model timestep; used here for STATS     [s]
         921             :
         922             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         923             :       threshold   ! Minimum value of x'^2                   [{x units}^2]
         924             :
         925             :     type (stats_metadata_type), intent(in) :: &
         926             :       stats_metadata
         927             :
         928             :     ! -------------------- InOut Variables --------------------
         929             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
         930             :       stats_zm
         931             :
         932             :     ! -------------------- Output Variable --------------------
         933             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
         934             :       xp2         ! Variance of x, x'^2 (momentum levels)   [{x units}^2]
         935             :
         936             :     ! -------------------- Local Variables --------------------
         937             :     integer :: i, k   ! Array index
         938             :
         939             :     integer :: &
         940             :       ixp2_cl
         941             :
         942             :     ! -------------------- Begin Code --------------------
         943             :
         944             :     !$acc data copyin( threshold ) &
         945             :     !$acc        copy( xp2 )
         946             :
         947    53610336 :     select case ( solve_type )
         948             :     case ( clip_wp2 )   ! wp2 clipping budget term
         949     8935056 :       ixp2_cl = stats_metadata%iwp2_cl
         950             :     case ( clip_rtp2 )   ! rtp2 clipping budget term
         951     8935056 :       ixp2_cl = stats_metadata%irtp2_cl
         952             :     case ( clip_thlp2 )   ! thlp2 clipping budget term
         953     8935056 :       ixp2_cl = stats_metadata%ithlp2_cl
         954             :     case ( clip_up2 )   ! up2 clipping budget term
         955     8935056 :       ixp2_cl = stats_metadata%iup2_cl
         956             :     case ( clip_vp2 )   ! vp2 clipping budget term
         957     8935056 :       ixp2_cl = stats_metadata%ivp2_cl
         958             :     case default   ! scalars are involved
         959    44675280 :       ixp2_cl = 0
         960             :     end select
         961             :
         962             :
         963    44675280 :     if ( stats_metadata%l_stats_samp ) then
         964             :       !$acc update host( xp2 )
         965           0 :       do i = 1, ngrdcol
         966           0 :         call stat_begin_update( nz, ixp2_cl, xp2(i,:) / dt, & ! intent(in)
         967           0 :                                 stats_zm(i) ) ! intent(inout)
         968             :       end do
         969             :     end if
         970             :
         971             :     ! Limit the value of x'^2 at threshold.
         972             :     ! The value of x'^2 at the surface (or lower boundary) is a set value that
         973             :     ! is determined elsewhere in a surface subroutine.  Thus, the variance
         974             :     ! clipping code does not need to be invoked at the lower boundary.
         975             :     ! Likewise, the value of x'^2 is set at the upper boundary, so the variance
         976             :     ! clipping code does not need to be invoked at the upper boundary.
         977             :     !
         978             :     ! charlass on 09/11/2013: I changed the clipping so that also the surface
         979             :     ! level is clipped. I did this because we discovered that there are slightly
         980             :     ! negative values in thlp2(1) and rtp2(1) when running quarter_ss case with
         981             :     ! WRF-CLUBB (see wrf:ticket:51#comment:33)
         982             :     !$acc parallel loop gang vector collapse(2) default(present)
         983  3797398800 :     do k = 1, nz-1, 1
         984 62706430800 :       do i = 1, ngrdcol
         985 62661755520 :         if ( xp2(i,k) < threshold(i,k) ) then
         986    60814029 :           xp2(i,k) = threshold(i,k)
         987             :         end if
         988             :       end do
         989             :     end do
         990             :     !$acc end parallel loop
         991             :
         992    44675280 :     if ( stats_metadata%l_stats_samp ) then
         993             :       !$acc update host( xp2 )
         994           0 :       do i = 1, ngrdcol
         995           0 :         call stat_end_update( nz, ixp2_cl, xp2(i,:) / dt, & ! intent(in)
         996           0 :                               stats_zm(i) ) ! intent(inout)
         997             :       end do
         998             :     end if
         999             :
        1000             :     !$acc end data
        1001             :
        1002    44675280 :     return
        1003             :
        1004             :   end subroutine clip_variance
        1005             :
        1006             :   !=============================================================================
        1007     8935056 :   subroutine clip_skewness( nz, ngrdcol, gr, dt, sfc_elevation, & ! intent(in)
        1008     8935056 :                             Skw_max_mag, wp2_zt,                & ! intent(in)
        1009             :                             l_use_wp3_lim_with_smth_Heaviside,  & ! intent(in)
        1010             :                             stats_metadata,                     & ! intent(in)
        1011     8935056 :                             stats_zt,                           & ! intent(inout)
        1012     8935056 :                             wp3 )                                 ! intent(out)
        1013             :
        1014             :     ! Description:
        1015             :     ! Clipping the value of w'^3 based on the skewness of w, Sk_w.
        1016             :     !
        1017             :     ! Aditionally, to prevent possible crashes due to wp3 growing too large,
        1018             :     ! abs(wp3) will be clipped to 100.
        1019             :     !
        1020             :     ! The skewness of w is:
        1021             :     !
        1022             :     ! Sk_w = w'^3 / (w'^2)^(3/2).
        1023             :     !
        1024             :     ! The value of Sk_w is limited to a range between an upper limit and a lower
        1025             :     ! limit.  The values of the limits depend on whether the level altitude is
        1026             :     ! within 100 meters of the surface.
        1027             :     !
        1028             :     ! For altitudes less than or equal to 100 meters above ground level (AGL):
        1029             :     !
        1030             :     ! -0.2_core_rknd*sqrt(2) <= Sk_w <= 0.2_core_rknd*sqrt(2);
        1031             :     !
        1032             :     ! while for all altitudes greater than 100 meters AGL:
        1033             :     !
        1034             :     ! -4.5_core_rknd <= Sk_w <= 4.5_core_rknd.
        1035             :     !
        1036             :     ! Therefore, there is an upper limit on w'^3, such that:
        1037             :     !
        1038             :     ! w'^3  <=  threshold_magnitude * (w'^2)^(3/2);
        1039             :     !
        1040             :     ! and a lower limit on w'^3, such that:
        1041             :     !
        1042             :     ! w'^3  >= -threshold_magnitude * (w'^2)^(3/2).
        1043             :     !
        1044             :     ! The values of w'^3 are found on the thermodynamic levels, while the values
        1045             :     ! of w'^2 are found on the momentum levels.  Therefore, the values of w'^2
        1046             :     ! are interpolated to the thermodynamic levels before being used to
        1047             :     ! calculate the upper and lower limits for w'^3.
        1048             :
        1049             :     ! References:
        1050             :     ! None
        1051             :     !-----------------------------------------------------------------------
        1052             :
        1053             :     use grid_class, only: &
        1054             :         grid ! Type
        1055             :
        1056             :     use clubb_precision, only: &
        1057             :         core_rknd ! Variable(s)
        1058             :
        1059             :     use stats_type_utilities, only: &
        1060             :         stat_begin_update,  & ! Procedure(s)
        1061             :         stat_end_update
        1062             :
        1063             :     use stats_variables, only: &
        1064             :         stats_metadata_type
        1065             :
        1066             :     use stats_type, only: stats ! Type
        1067             :
        1068             :     implicit none
        1069             :
        1070             :     ! ----------------------- Input Variables -----------------------
        1071             :     integer, intent(in) :: &
        1072             :       nz, &
        1073             :       ngrdcol
        1074             :
        1075             :     type (grid), target, intent(in) :: gr
        1076             :
        1077             :     real( kind = core_rknd ), intent(in) :: &
        1078             :       dt               ! Model timestep; used here for STATS        [s]
        1079             :
        1080             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
        1081             :       sfc_elevation,    & ! Elevation of ground level                  [m AMSL]
        1082             :       Skw_max_mag         ! Maximum allowable magnitude of Skewness    [-]
        1083             :
        1084             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1085             :       wp2_zt           ! w'^2 interpolated to thermodyamic levels   [m^2/s^2]
        1086             :
        1087             :     ! Flag to activate modifications on wp3 limiters for convergence test
        1088             :     ! (use smooth Heaviside 'Preskin' function in the calculation of
        1089             :     ! clip_skewness for wp3)
        1090             :     logical, intent(in):: &
        1091             :       l_use_wp3_lim_with_smth_Heaviside
        1092             :
        1093             :     type (stats_metadata_type), intent(in) :: &
        1094             :       stats_metadata
        1095             :
        1096             :     ! ----------------------- Input/Output Variables -----------------------
        1097             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
        1098             :       stats_zt
        1099             :
        1100             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
        1101             :       wp3              ! w'^3 (thermodynamic levels)                [m^3/s^3]
        1102             :
        1103             :     ! ----------------------- Local Variables -----------------------
        1104             :     integer :: i
        1105             :
        1106             :     ! ----------------------- Begin Code -----------------------
        1107             :
        1108             :     !$acc data copyin( gr, gr%zt, &
        1109             :     !$acc              sfc_elevation, wp2_zt ) &
        1110             :     !$acc        copy( wp3 )
        1111             :
        1112     8935056 :     if ( stats_metadata%l_stats_samp ) then
        1113             :
        1114             :       !$acc update host( wp3 )
        1115             :
        1116           0 :       do i = 1, ngrdcol
        1117           0 :         call stat_begin_update( nz, stats_metadata%iwp3_cl, wp3(i,:) / dt, & ! intent(in)
        1118           0 :                                 stats_zt(i) ) ! intent(inout)
        1119             :       end do
        1120             :     end if
        1121             :
        1122             :     call clip_skewness_core( nz, ngrdcol, gr, sfc_elevation,    & ! intent(in)
        1123             :                              Skw_max_mag, wp2_zt,               & ! intent(in)
        1124             :                              l_use_wp3_lim_with_smth_Heaviside, & ! intent(in)
        1125     8935056 :                              wp3 )                                ! intent(inout)
        1126             :
        1127     8935056 :     if ( stats_metadata%l_stats_samp ) then
        1128             :
        1129             :       !$acc update host( wp3 )
        1130             :
        1131           0 :       do i = 1, ngrdcol
        1132           0 :         call stat_end_update( nz, stats_metadata%iwp3_cl, wp3(i,:) / dt, & ! intent(in)
        1133           0 :                               stats_zt(i) ) ! intent(inout)
        1134             :       end do
        1135             :     end if
        1136             :
        1137             :     !$acc end data
        1138             :
        1139     8935056 :     return
        1140             :
        1141             :   end subroutine clip_skewness
        1142             :
        1143             : !=============================================================================
        1144     8935056 :   subroutine clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, &
        1145     8935056 :                                  Skw_max_mag, wp2_zt, &
        1146             :                                  l_use_wp3_lim_with_smth_Heaviside, &
        1147     8935056 :                                  wp3 )
        1148             :
        1149             :     use grid_class, only: &
        1150             :         grid ! Type
        1151             :
        1152             :     use clubb_precision, only: &
        1153             :         core_rknd ! Variable(s)
        1154             :
        1155             :     use advance_helper_module, only: &
        1156             :         smooth_heaviside_peskin
        1157             :
        1158             :     implicit none
        1159             :
        1160             :     !----------------------- Input Variables -----------------------
        1161             :     integer, intent(in) :: &
        1162             :       nz, &
        1163             :       ngrdcol
        1164             :
        1165             :     type (grid), target, intent(in) :: gr
        1166             :
        1167             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
        1168             :       sfc_elevation,  & ! Elevation of ground level                  [m AMSL]
        1169             :       Skw_max_mag       ! Maximum allowable magnitude of Skewness    [-]
        1170             :
        1171             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1172             :       wp2_zt           ! w'^2 interpolated to thermodyamic levels   [m^2/s^2]
        1173             :
        1174             :     ! Flag to activate modifications on wp3 limiters for convergence test
        1175             :     ! (use smooth Heaviside 'Preskin' function in the calculation of clip_skewness for wp3)
        1176             :     logical, intent(in):: &
        1177             :       l_use_wp3_lim_with_smth_Heaviside
        1178             :
        1179             :     !----------------------- Input/Output Variables -----------------------
        1180             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
        1181             :       wp3              ! w'^3 (thermodynamic levels)                [m^3/s^3]
        1182             :
        1183             :     !----------------------- Local Variables -----------------------
        1184             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
        1185    17870112 :       wp2_zt_cubed, & ! Variance of vertical velocity cubed (w^2_{zt}^3)   [m^6/s^6]
        1186    17870112 :       wp3_lim_sqd     ! Keeps absolute value of Sk_w from becoming > limit [m^6/s^6]
        1187             :
        1188             :     integer :: i, k       ! Vertical array index.
        1189             :
        1190             :     real( kind = core_rknd ), parameter :: &
        1191             :       wp3_max = 100._core_rknd ! Threshold for wp3 [m^3/s^3]
        1192             :
        1193             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
        1194    17870112 :       zagl_thresh, & ! temporatory array
        1195    17870112 :       H_zagl ! Heaviside function for clippings of wp3_lim_sqd
        1196             :
        1197             :     !----------------------- Begin Code-----------------------
        1198             :
        1199             :     !$acc enter data create( wp2_zt_cubed, wp3_lim_sqd, zagl_thresh, H_zagl )
        1200             :
        1201             :     ! Compute the upper and lower limits of w'^3 at every level,
        1202             :     ! based on the skewness of w, Sk_w, such that:
        1203             :     ! Sk_w = w'^3 / (w'^2)^(3/2);
        1204             :     ! -4.5 <= Sk_w <= 4.5;
        1205             :     ! or, if the level altitude is within 100 meters of the surface,
        1206             :     ! -0.2*sqrt(2) <= Sk_w <= 0.2*sqrt(2).
        1207             :
        1208             :     ! The normal magnitude limit of skewness of w in the CLUBB code is 4.5.
        1209             :     ! However, according to Andre et al. (1976b & 1978), wp3 should not exceed
        1210             :     ! [2*(wp2^3)]^(1/2) at any level.  However, this term should be multiplied
        1211             :     ! by 0.2 close to the surface to include surface effects.  There already is
        1212             :     ! a wp3 clipping term in place for all other altitudes, but this term will
        1213             :     ! be included for the surface layer only.  Therefore, the lowest level wp3
        1214             :     ! should not exceed 0.2 * sqrt(2) * wp2^(3/2).  Brian Griffin.  12/18/05.
        1215             :
        1216             :     ! To lower compute time, we squared both sides of the equation and compute
        1217             :     ! wp2^3 only once. -dschanen 9 Oct 2008
        1218             :     !$acc parallel loop gang vector collapse(2) default(present)
        1219   768414816 :     do k = 1, nz
        1220 12690480816 :       do i = 1, ngrdcol
        1221 12681545760 :         wp2_zt_cubed(i,k) = wp2_zt(i,k)**3
        1222             :       end do
        1223             :     end do
        1224             :     !$acc end parallel loop
        1225             :
        1226     8935056 :     if ( l_use_wp3_lim_with_smth_Heaviside ) then
        1227             :
        1228             :       !implement a smoothed Heaviside function to avoid discontinuities
        1229             :       !$acc parallel loop gang vector collapse(2) default(present)
        1230           0 :       do k = 1, nz
        1231           0 :         do i = 1, ngrdcol
        1232           0 :           zagl_thresh(i,k) = ( gr%zt(i,k) - sfc_elevation(i) ) /  100.0_core_rknd
        1233           0 :           zagl_thresh(i,k) = zagl_thresh(i,k)  - 1.0_core_rknd
        1234             :         end do
        1235             :       end do
        1236             :       !$acc end parallel loop
        1237             :
        1238           0 :       H_zagl(:,:) = smooth_heaviside_peskin(nz, ngrdcol, zagl_thresh(:,:), 0.6_core_rknd)
        1239             :
        1240             :       !$acc parallel loop gang vector collapse(2) default(present)
        1241           0 :       do k = 1, nz
        1242           0 :         do i = 1, ngrdcol
        1243           0 :           wp3_lim_sqd(i,k) = wp2_zt_cubed(i,k)   &
        1244             :                               * ( H_zagl(i,k) * Skw_max_mag(i)**2   &
        1245             :                                   + (1.0_core_rknd - H_zagl(i,k)) &
        1246           0 :                                      * 0.0021_core_rknd *Skw_max_mag(i)**2 )
        1247             :         end do
        1248             :       end do
        1249             :      !$acc end parallel loop
        1250             :
        1251             :     else ! default method
        1252             :
        1253             :       !$acc parallel loop gang vector collapse(2) default(present)
        1254   768414816 :       do k = 1, nz
        1255 12690480816 :         do i = 1, ngrdcol
        1256 12681545760 :           if ( gr%zt(i,k) - sfc_elevation(i) <= 100.0_core_rknd ) then ! Clip for 100 m. AGL.
        1257             :            !wp3_upper_lim(k) =  0.2_core_rknd * sqrt_2 * wp2_zt(k)**(3.0_core_rknd/2.0_core_rknd)
        1258             :            !wp3_lower_lim(k) = -0.2_core_rknd * sqrt_2 * wp2_zt(k)**(3.0_core_rknd/2.0_core_rknd)
        1259   420778800 :             wp3_lim_sqd(i,k) = 0.0021_core_rknd * Skw_max_mag(i)**2 * wp2_zt_cubed(i,k)
        1260             :           else                          ! Clip skewness consistently with a.
        1261             :            !wp3_upper_lim(k) =  4.5_core_rknd * wp2_zt(k)**(3.0_core_rknd/2.0_core_rknd)
        1262             :            !wp3_lower_lim(k) = -4.5_core_rknd * wp2_zt(k)**(3.0_core_rknd/2.0_core_rknd)
        1263 11501287200 :             wp3_lim_sqd(i,k) = Skw_max_mag(i)**2 * wp2_zt_cubed(i,k) ! Skw_max_mag = 4.5_core_rknd^2
        1264             :           endif
        1265             :         end do
        1266             :       end do
        1267             :       !$acc end parallel loop
        1268             :
        1269             :     end if
        1270             :
        1271             :     ! Clipping for w'^3 at an upper and lower limit corresponding with
        1272             :     ! the appropriate value of Sk_w.
        1273             :     !$acc parallel loop gang vector collapse(2) default(present)
        1274   768414816 :     do k = 1, nz
        1275 12690480816 :       do i = 1, ngrdcol
        1276             :         ! Set the magnitude to the wp3 limit and apply the sign of the current wp3
        1277 12681545760 :         if ( wp3(i,k)**2 > wp3_lim_sqd(i,k) ) then
        1278   206449814 :           wp3(i,k) = sign( sqrt( wp3_lim_sqd(i,k) ), wp3(i,k) )
        1279             :         end if
        1280             :       end do
        1281             :     end do
        1282             :     !$acc end parallel loop
        1283             :
        1284             :     ! Clipping abs(wp3) to 100. This keeps wp3 from growing too large in some
        1285             :     ! deep convective cases, which helps prevent these cases from blowing up.
        1286             :     !$acc parallel loop gang vector collapse(2) default(present)
        1287   768414816 :     do k = 1, nz
        1288 12690480816 :       do i = 1, ngrdcol
        1289 12681545760 :         if ( abs(wp3(i,k)) > wp3_max ) then
        1290          32 :           wp3(i,k) = sign( wp3_max, wp3(i,k) ) ! Known magic number
        1291             :         end if
        1292             :       end do
        1293             :     end do
        1294             :     !$acc end parallel loop
        1295             :
        1296             :     !$acc exit data delete( wp2_zt_cubed, wp3_lim_sqd, zagl_thresh, H_zagl )
        1297             :
        1298     8935056 :   end subroutine clip_skewness_core
        1299             :
        1300             : !===============================================================================
        1301             :
        1302             : end module clip_explicit

Generated by: LCOV version 1.14