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 17:57:11 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    34070823 :           xpyp_chnge(i,k) = xpyp_bound - xpyp(i,k)
     619    34070823 :           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    30885521 :           xpyp_chnge(i,k) = -xpyp_bound - xpyp(i,k)
     626    30885521 :           xpyp(i,k) = -xpyp_bound 
     627             : 
     628             :         else
     629             : 
     630 >15127*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    57712676 :           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   201653678 :           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           4 :           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