LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - sfc_varnce_module.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 68 | 140 | 48.6 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 1 | 100.0 %

              Line data    Source code
           1             : !-------------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module sfc_varnce_module
           5             :
           6             :   implicit none
           7             :
           8             :   private ! Default to private
           9             :
          10             :   public :: calc_sfc_varnce
          11             :
          12             :   contains
          13             :
          14             :   !=============================================================================
          15     8935056 :   subroutine calc_sfc_varnce( nz, ngrdcol, sclr_dim, sclr_idx, &
          16     8935056 :                               gr, dt, sfc_elevation, &
          17     8935056 :                               upwp_sfc, vpwp_sfc, wpthlp, wprtp_sfc, &
          18     8935056 :                               um, vm, Lscale_up, wpsclrp_sfc, &
          19     8935056 :                               lhs_splat_wp2, tau_zm, &
          20             :                               !wp2_splat_sfc, tau_zm_sfc, &
          21             :                               l_vary_convect_depth, &
          22     8935056 :                               up2_sfc_coef, &
          23     8935056 :                               a_const, &
          24             :                               stats_metadata, &
          25     8935056 :                               stats_zm, &
          26     8935056 :                               wp2, up2, vp2, &
          27     8935056 :                               thlp2, rtp2, rtpthlp, &
          28     8935056 :                               sclrp2, sclrprtp,  sclrpthlp )
          29             :
          30             :     ! Description:
          31             :     ! This subroutine computes estimate of the surface thermodynamic and wind
          32             :     ! component second order moments.
          33             :
          34             :     ! References:
          35             :     ! Andre, J. C., G. De Moor, P. Lacarrere, G. Therry, and R. Du Vachat, 1978.
          36             :     !   Modeling the 24-Hour Evolution of the Mean and Turbulent Structures of
          37             :     !   the Planetary Boundary Layer.  J. Atmos. Sci., 35, 1861 -- 1883.
          38             :     !-----------------------------------------------------------------------
          39             :
          40             :     use grid_class, only: &
          41             :       grid
          42             :
          43             :     use parameters_model, only:  &
          44             :         T0 ! Variable(s)
          45             :
          46             :     use constants_clubb, only: &
          47             :         four,       & ! Variable(s)
          48             :         two,        &
          49             :         one,        &
          50             :         two_thirds, &
          51             :         one_third,  &
          52             :         one_fourth, &
          53             :         zero,       &
          54             :         grav,       &
          55             :         eps,        &
          56             :         w_tol_sqd,  &
          57             :         thl_tol,    &
          58             :         rt_tol,     &
          59             :         max_mag_correlation_flux, &
          60             :         fstderr,    &
          61             :         wp2_max
          62             :
          63             :     use numerical_check, only: &
          64             :         sfc_varnce_check ! Procedure
          65             :
          66             :     use error_code, only: &
          67             :         clubb_at_least_debug_level,  & ! Procedure
          68             :         err_code,                    & ! Error Indicator
          69             :         clubb_fatal_error              ! Constant
          70             :
          71             :     use array_index, only: &
          72             :         sclr_idx_type
          73             :
          74             :     use stats_type, only: &
          75             :         stats ! Type
          76             :
          77             :     use stats_type_utilities, only: &
          78             :         stat_end_update_pt,   & ! Procedure(s)
          79             :         stat_begin_update_pt, &
          80             :         stat_update_var_pt
          81             :
          82             :     use stats_variables, only: &
          83             :         stats_metadata_type
          84             :
          85             :     use clubb_precision, only: &
          86             :         core_rknd ! Variable(s)
          87             :
          88             :     implicit none
          89             :
          90             :     ! Constant Parameters
          91             :
          92             :     ! Logical for Andre et al., 1978 parameterization.
          93             :     logical, parameter :: l_andre_1978 = .false.
          94             :
          95             :     real( kind = core_rknd ), parameter ::  &
          96             :       z_const = one, & ! Defined height of 1 meter                [m]
          97             :         ! Vince Larson increased ufmin to stabilize arm_97.  24 Jul 2007
          98             :         ! ufmin = 0.0001_core_rknd, &
          99             :       ufmin = 0.01_core_rknd, & ! Minimum allowable value of u*   [m/s]
         100             :         ! End Vince Larson's change.
         101             :         ! Vince Larson changed in order to make correlations between [-1,1].  31 Jan 2008.
         102             :         ! sclr_var_coef = 0.25_core_rknd, & ! This value is made up! - Vince Larson 12 Jul 2005
         103             :       sclr_var_coef = 0.4_core_rknd,  & ! This value is made up! - Vince Larson 12 Jul 2005
         104             :         ! End Vince Larson's change
         105             :         ! Vince Larson reduced surface spike in scalar variances associated
         106             :         ! w/ Andre et al. 1978 scheme
         107             :       reduce_coef   = 0.2_core_rknd
         108             :
         109             :     !-------------------------- Input Variables --------------------------
         110             :     integer, intent(in) :: &
         111             :       nz, &
         112             :       ngrdcol, &
         113             :       sclr_dim
         114             :
         115             :     type (sclr_idx_type), intent(in) :: &
         116             :       sclr_idx
         117             :
         118             :     type(grid), intent(in) :: &
         119             :       gr
         120             :
         121             :     real( kind = core_rknd ) ::  &
         122             :       dt
         123             :
         124             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
         125             :       sfc_elevation, &
         126             :       upwp_sfc,         & ! Surface u momentum flux, <u'w'>|_sfc   [m^2/s^2]
         127             :       vpwp_sfc,         & ! Surface v momentum flux, <v'w'>|_sfc   [m^2/s^2]
         128             :       wprtp_sfc           ! Surface moisture flux, <w'rt'>|_sfc    [kg/kg m/s]
         129             :
         130             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         131             :       wpthlp,        & ! Surface thetal flux, <w'thl'>|_sfc     [K m/s]
         132             :       um,            & ! Surface u wind component, <u>          [m/s]
         133             :       vm,            & ! Surface v wind component, <v>          [m/s]
         134             :       Lscale_up,     & ! Upward component of Lscale at surface  [m]
         135             :       lhs_splat_wp2, &
         136             :       !wp2_splat,    & ! Tendency of <w'^2> due to splatting of eddies at zm(1) [m^2/s^3]
         137             :       tau_zm           ! Turbulent dissipation time at level zm(1)  [s]
         138             :
         139             :
         140             :     real( kind = core_rknd ), dimension(ngrdcol,sclr_dim), intent(in) ::  &
         141             :       wpsclrp_sfc    ! Passive scalar flux, <w'sclr'>|_sfc   [units m/s]
         142             :
         143             :     logical, intent(in) :: &
         144             :       l_vary_convect_depth
         145             :
         146             :     real( kind = core_rknd ), dimension(ngrdcol) :: &
         147             :       up2_sfc_coef,   & ! CLUBB tunable parameter up2_sfc_coef   [-]
         148             :       a_const           ! Coefficient in front of wp2, up2, and vp2
         149             :
         150             :     type (stats_metadata_type), intent(in) :: &
         151             :       stats_metadata
         152             :
         153             :     !-------------------------- InOut Variables --------------------------
         154             :     type (stats), target, intent(inout), dimension(ngrdcol) :: &
         155             :       stats_zm
         156             :
         157             :     !-------------------------- Output Variables --------------------------
         158             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) ::  &
         159             :       wp2,     & ! Surface variance of w, <w'^2>|_sfc            [m^2/s^2]
         160             :       up2,     & ! Surface variance of u, <u'^2>|_sfc            [m^2/s^2]
         161             :       vp2,     & ! Surface variance of v, <v'^2>|_sfc            [m^2/s^2]
         162             :       thlp2,   & ! Surface variance of theta-l, <thl'^2>|_sfc    [K^2]
         163             :       rtp2,    & ! Surface variance of rt, <rt'^2>|_sfc          [(kg/kg)^2]
         164             :       rtpthlp    ! Surface covariance of rt and theta-l          [kg K/kg]
         165             :
         166             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(inout) ::  &
         167             :       sclrp2,    & ! Surface variance of passive scalar            [units^2]
         168             :       sclrprtp,  & ! Surface covariance of pssv scalar and rt  [units kg/kg]
         169             :       sclrpthlp    ! Surface covariance of pssv scalar and theta-l [units K]
         170             :
         171             :     !-------------------------- Local Variables --------------------------
         172             :     real( kind = core_rknd ), dimension(ngrdcol) ::  &
         173    17870112 :       uf,               & ! Surface friction vel., u*, in older version   [m/s]
         174    17870112 :       depth_pos_wpthlp, & ! Thickness of the layer near the surface with wpthlp > 0 [m]
         175    17870112 :       min_wp2_sfc_val     ! Minimum value of wp2_sfc that guarantees  [m^2/s^2]
         176             :                           ! correlation of (w,rt) and (w,thl) is within (-1,1)
         177             :
         178             :     ! Variables for Andre et al., 1978 parameterization.
         179             :     real( kind = core_rknd ), dimension(ngrdcol) :: &
         180     8935056 :       um_sfc_sqd, & ! Surface value of <u>^2                           [m^2/s^2]
         181     8935056 :       vm_sfc_sqd, & ! Surface value of <v>^2                           [m^2/s^2]
         182     8935056 :       usp2_sfc,   & ! u_s (vector oriented w/ mean sfc. wind) variance [m^2/s^2]
         183     8935056 :       vsp2_sfc,   & ! v_s (vector perpen. to mean sfc. wind) variance  [m^2/s^2]
         184     8935056 :       ustar,      & ! Surface friction velocity, u*                             [m/s]
         185     8935056 :       zeta,       & ! Dimensionless height z_const/Lngth, where z_const = 1 m.  [-]
         186    17870112 :       wp2_splat_sfc_correction  ! Reduction in wp2_sfc due to splatting [m^2/s^2]
         187             :
         188             :     real( kind = core_rknd ) ::  &
         189             :       ustar2, & ! Square of surface friction velocity, u*       [m^2/s^2]
         190             :       wstar     ! Convective velocity, w*                       [m/s]
         191             :
         192             :     real( kind = core_rknd ) :: &
         193             :       Lngth    ! Monin-Obukhov length [m]
         194             :
         195             :     integer :: i, k, sclr ! Loop index
         196             :
         197             :     !-------------------------- Begin Code --------------------------
         198             :
         199             :     !$acc enter data create( uf, depth_pos_wpthlp, min_wp2_sfc_val, &
         200             :     !$acc                    um_sfc_sqd, vm_sfc_sqd, usp2_sfc, vsp2_sfc, &
         201             :     !$acc                    ustar, zeta, wp2_splat_sfc_correction )
         202             :
         203             :     ! Reflect surface varnce changes in budget
         204     8935056 :     if ( stats_metadata%l_stats_samp ) then
         205             :
         206             :       !$acc update host( wp2, up2, vp2, thlp2, rtp2, rtpthlp )
         207             :
         208           0 :       do i = 1, ngrdcol
         209             :         call stat_begin_update_pt( stats_metadata%ithlp2_sf, 1,      & ! intent(in)
         210           0 :                                    thlp2(i,1) / dt,   & ! intent(in)
         211           0 :                                    stats_zm(i) )        ! intent(inout)
         212             :
         213             :         call stat_begin_update_pt( stats_metadata%irtp2_sf, 1,       & ! intent(in)
         214             :                                    rtp2(i,1) / dt,    & ! intent(in)
         215           0 :                                    stats_zm(i) )        ! intent(inout)
         216             :
         217             :         call stat_begin_update_pt( stats_metadata%irtpthlp_sf, 1,    & ! intent(in)
         218             :                                    rtpthlp(i,1) / dt, & ! intent(in)
         219           0 :                                    stats_zm(i) )        ! intent(inout)
         220             :
         221             :         call stat_begin_update_pt( stats_metadata%iup2_sf, 1,        & ! intent(in)
         222             :                                    up2(i,1) / dt,     & ! intent(in)
         223           0 :                                    stats_zm(i) )        ! intent(inout)
         224             :
         225             :         call stat_begin_update_pt( stats_metadata%ivp2_sf, 1,        & ! intent(in)
         226             :                                    vp2(i,1) / dt,     & ! intent(in)
         227           0 :                                    stats_zm(i) ) ! intent(inout)
         228             :
         229             :         call stat_begin_update_pt( stats_metadata%iwp2_sf, 1,        & ! intent(in)
         230             :                                    wp2(i,1) / dt,     & ! intent(in)
         231           0 :                                    stats_zm(i) )        ! intent(inout)
         232             :       end do
         233             :     end if
         234             :
         235             :     !$acc parallel loop gang vector default(present)
         236   149194656 :     do i = 1, ngrdcol
         237             :       ! Find thickness of layer near surface with positive heat flux.
         238             :       ! This is used when l_vary_convect_depth=.true. in order to determine wp2.
         239   149194656 :       if ( wpthlp(i,1) <= zero ) then
         240    46511616 :          depth_pos_wpthlp(i) = one ! When sfc heat flux is negative, set depth to 1 m.
         241             :       else ! When sfc heat flux is positive, march up sounding until wpthlp 1st becomes negative.
         242             :          k = 1
         243   571111866 :          do while ( wpthlp(i,k) > zero .and. (gr%zm(i,k)-sfc_elevation(i)) < 1000._core_rknd )
         244   477363882 :             k = k + 1
         245             :         end do
         246    93747984 :         depth_pos_wpthlp(i) = max( one, gr%zm(i,k)-sfc_elevation(i) )
         247             :       end if
         248             :     end do
         249             :     !$acc end parallel loop
         250             :
         251             :     ! a_const used to be set here; now it is a tunable parameter
         252             :     !if ( .not. l_vary_convect_depth ) then
         253             :     !   a_const = 1.8_core_rknd
         254             :     !else
         255             :     !   a_const = 0.6_core_rknd
         256             :     !end if
         257             :
         258             :     if ( l_andre_1978 ) then
         259             :
         260             :       ! Calculate <u>^2 and <v>^2.
         261             :       !$acc parallel loop gang vector default(present)
         262             :       do i = 1, ngrdcol
         263             :         um_sfc_sqd(i) = um(i,2)**2
         264             :         vm_sfc_sqd(i) = vm(i,2)**2
         265             :       end do
         266             :       !$acc end parallel loop
         267             :
         268             :       ! Calculate surface friction velocity, u*.
         269             :       !$acc parallel loop gang vector default(present)
         270             :       do i = 1, ngrdcol
         271             :         ustar(i) = max( ( upwp_sfc(i)**2 + vpwp_sfc(i)**2 )**(one_fourth), ufmin )
         272             :       end do
         273             :       !$acc end parallel loop
         274             :
         275             :       !$acc parallel loop gang vector default(present)
         276             :       do i = 1, ngrdcol
         277             :         if ( abs(wpthlp(i,1)) > eps) then
         278             :
         279             :           ! Find Monin-Obukhov Length (Andre et al., 1978, p. 1866).
         280             :           Lngth = - ( ustar(i)**3 ) &
         281             :                     / ( 0.35_core_rknd * (one/T0) * grav * wpthlp(i,1) )
         282             :
         283             :           ! Find the value of dimensionless height zeta
         284             :           ! (Andre et al., 1978, p. 1866).
         285             :           zeta(i) = z_const / Lngth
         286             :
         287             :        else ! wpthlp = 0
         288             :
         289             :           ! The value of Monin-Obukhov length is +inf when ustar < 0 and -inf
         290             :           ! when ustar > 0.  Either way, zeta = 0.
         291             :           zeta(i) = zero
         292             :
         293             :         endif ! wpthlp /= 0
         294             :       end do
         295             :       !$acc end parallel loop
         296             :
         297             :       ! Andre et al, 1978, Eq. 29.
         298             :       ! Notes:  1) "reduce_coef" is a reduction coefficient intended to make
         299             :       !            the values of rtp2, thlp2, and rtpthlp smaller at the
         300             :       !            surface.
         301             :       !         2) With the reduction coefficient having a value of 0.2, the
         302             :       !            surface correlations of both w & rt and w & thl have a value
         303             :       !            of about 0.845.  These correlations are greater if zeta < 0.
         304             :       !            The correlations have a value greater than 1 if
         305             :       !            zeta <= -0.212.
         306             :       !         3) The surface correlation of rt & thl is 1.
         307             :       ! Brian Griffin; February 2, 2008.
         308             :       !$acc parallel loop gang vector default(present)
         309             :       do i = 1, ngrdcol
         310             :         if ( zeta(i) < zero ) then
         311             :
         312             :           thlp2(i,1)   = reduce_coef  &
         313             :                         * ( wpthlp(i,1)**2 / ustar(i)**2 ) &
         314             :                         * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
         315             :
         316             :           rtp2(i,1)    = reduce_coef  &
         317             :                         * ( wprtp_sfc(i)**2 / ustar(i)**2 ) &
         318             :                         * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
         319             :
         320             :           rtpthlp(i,1) = reduce_coef  &
         321             :                         * ( wprtp_sfc(i) * wpthlp(i,1) / ustar(i)**2 ) &
         322             :                         * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
         323             :
         324             :           wp2(i,1)    = ( ustar(i)**2 ) &
         325             :                           * ( 1.75_core_rknd + two * (-zeta(i))**(two_thirds) )
         326             :
         327             :        else
         328             :
         329             :           thlp2(i,1)   = reduce_coef  &
         330             :                         * four * ( wpthlp(i,1)**2 / ustar(i)**2 )
         331             :
         332             :           rtp2(i,1)    = reduce_coef  &
         333             :                         * four * ( wprtp_sfc(i)**2 / ustar(i)**2 )
         334             :
         335             :           rtpthlp(i,1) = reduce_coef  &
         336             :                         * four * ( wprtp_sfc(i) * wpthlp(i,1) / ustar(i)**2 )
         337             :
         338             :           wp2(i,1)     = 1.75_core_rknd * ustar(i)**2
         339             :
         340             :         endif
         341             :       end do
         342             :       !$acc end parallel loop
         343             :
         344             :       !$acc parallel loop gang vector default(present)
         345             :       do i = 1, ngrdcol
         346             :         thlp2(i,1) = max( thl_tol**2, thlp2(i,1) )
         347             :         rtp2(i,1) = max( rt_tol**2, rtp2(i,1) )
         348             :       end do
         349             :       !$acc end parallel loop
         350             :
         351             :       ! Calculate wstar following Andre et al., 1978, p. 1866.
         352             :       ! w* = ( ( 1 / T0 ) * g * <w'thl'>|_sfc * z_i )^(1/3);
         353             :       ! where z_i is the height of the mixed layer.  The value of CLUBB's
         354             :       ! upward component of mixing length, Lscale_up, at the surface will be
         355             :       ! used as z_i.
         356             :       !$acc parallel loop gang vector default(present)
         357             :       do i = 1, ngrdcol
         358             :         wstar = ( (one/T0) * grav * wpthlp(i,1) * Lscale_up(i,2) )**(one_third)
         359             :
         360             :         ! Andre et al., 1978, Eq. 29.
         361             :         ! Andre et al. (1978) defines horizontal wind surface variances in terms
         362             :         ! of orientation with the mean surface wind.  The vector u_s is the wind
         363             :         ! vector oriented with the mean surface wind.  The vector v_s is the wind
         364             :         ! vector oriented perpendicular to the mean surface wind.  Thus, <u_s> is
         365             :         ! equal to the mean surface wind (both in speed and direction), and <v_s>
         366             :         ! is 0.  Equation 29 gives the formula for the variance of u_s, which is
         367             :         ! <u_s'^2> (usp2_sfc in the code), and the formula for the variance of
         368             :         ! v_s, which is <v_s'^2> (vsp2_sfc in the code).
         369             :         if ( wpthlp(i,1) > zero ) then
         370             :
         371             :           usp2_sfc(i) = four * ustar(i)**2 + 0.3_core_rknd * wstar**2
         372             :
         373             :           vsp2_sfc(i) = 1.75_core_rknd * ustar(i)**2 + 0.3_core_rknd * wstar**2
         374             :
         375             :         else
         376             :
         377             :           usp2_sfc(i) = four * ustar(i)**2
         378             :
         379             :           vsp2_sfc(i) = 1.75_core_rknd * ustar(i)**2
         380             :
         381             :         endif
         382             :       end do
         383             :       !$acc end parallel loop
         384             :
         385             :       ! Add effect of vertical compression of eddies on horizontal gustiness.
         386             :       ! First, ensure that wp2 does not make the correlation
         387             :       !   of (w,thl) or (w,rt) outside (-1,1).
         388             :       ! Perhaps in the future we should also ensure that the correlations
         389             :       !   of (w,u) and (w,v) are not outside (-1,1).
         390             :       !$acc parallel loop gang vector default(present)
         391             :       do i = 1, ngrdcol
         392             :         min_wp2_sfc_val(i) &
         393             :                = max( w_tol_sqd, &
         394             :                       wprtp_sfc(i)**2 / ( rtp2(i,1) * max_mag_correlation_flux**2 ), &
         395             :                       wpthlp(i,1)**2 / ( thlp2(i,1) * max_mag_correlation_flux**2 ) )
         396             :       end do
         397             :       !$acc end parallel loop
         398             :
         399             :       !$acc parallel loop gang vector default(present)
         400             :       do i = 1, ngrdcol
         401             :         if ( wp2(i,1) - tau_zm(i,1) * lhs_splat_wp2(i,1) * wp2(i,1) &
         402             :              < min_wp2_sfc_val(i) ) then
         403             :         !if ( wp2(i,1) + tau_zm(i,1) * wp2_splat(i,1) < min_wp2_sfc_val(i) ) then
         404             :                               ! splatting correction drives wp2 to overly small value
         405             :           wp2_splat_sfc_correction(i) = -wp2(i,1) + min_wp2_sfc_val(i)
         406             :           wp2(i,1) = min_wp2_sfc_val(i)
         407             :         else
         408             :           wp2_splat_sfc_correction(i) = tau_zm(i,1) * lhs_splat_wp2(i,1) * wp2(i,1)
         409             :           !wp2_splat_sfc_correction(i) = tau_zm(i,1) * wp2_splat(i,1)
         410             :           wp2(i,1) = wp2(i,1) + wp2_splat_sfc_correction(i)
         411             :         end if
         412             :       end do
         413             :       !$acc end parallel loop
         414             :
         415             :       !$acc parallel loop gang vector default(present)
         416             :       do i = 1, ngrdcol
         417             :         usp2_sfc(i) = usp2_sfc(i) - 0.5_core_rknd * wp2_splat_sfc_correction(i)
         418             :         vsp2_sfc(i) = vsp2_sfc(i) - 0.5_core_rknd * wp2_splat_sfc_correction(i)
         419             :       end do
         420             :       !$acc end parallel loop
         421             :
         422             :       ! Variance of u, <u'^2>, at the surface can be found from <u_s'^2>,
         423             :       ! <v_s'^2>, and mean winds (at the surface) <u> and <v>, such that:
         424             :       !    <u'^2>|_sfc = <u_s'^2> * [ <u>^2 / ( <u>^2 + <v>^2 ) ]
         425             :       !                  + <v_s'^2> * [ <v>^2 / ( <u>^2 + <v>^2 ) ];
         426             :       ! where <u>^2 + <v>^2 /= 0.
         427             :       !$acc parallel loop gang vector default(present)
         428             :       do i = 1, ngrdcol
         429             :         up2(i,1) = usp2_sfc(i) * ( um_sfc_sqd(i) / max( um_sfc_sqd(i) + vm_sfc_sqd(i) , eps ) )  &
         430             :                    + vsp2_sfc(i) * ( vm_sfc_sqd(i) / max( um_sfc_sqd(i) + vm_sfc_sqd(i) , eps ) )
         431             :       end do
         432             :       !$acc end parallel loop
         433             :
         434             :       ! Variance of v, <v'^2>, at the surface can be found from <u_s'^2>,
         435             :       ! <v_s'^2>, and mean winds (at the surface) <u> and <v>, such that:
         436             :       !    <v'^2>|_sfc = <v_s'^2> * [ <u>^2 / ( <u>^2 + <v>^2 ) ]
         437             :       !                  + <u_s'^2> * [ <v>^2 / ( <u>^2 + <v>^2 ) ];
         438             :       ! where <u>^2 + <v>^2 /= 0.
         439             :       !$acc parallel loop gang vector default(present)
         440             :       do i = 1, ngrdcol
         441             :         vp2(i,1) = vsp2_sfc(i) * ( um_sfc_sqd(i) / max( um_sfc_sqd(i) + vm_sfc_sqd(i) , eps ) )  &
         442             :                    + usp2_sfc(i) * ( vm_sfc_sqd(i) / max( um_sfc_sqd(i) + vm_sfc_sqd(i) , eps ) )
         443             :       end do
         444             :       !$acc end parallel loop
         445             :
         446             :       ! Passive scalars
         447             :       if ( sclr_dim > 0 ) then
         448             :
         449             :         !$acc parallel loop gang vector collapse(2) default(present)
         450             :         do sclr = 1, sclr_dim
         451             :           do i = 1, ngrdcol
         452             :
         453             :             ! Notes:  1) "reduce_coef" is a reduction coefficient intended to
         454             :             !            make the values of sclrprtp, sclrpthlp, and sclrp2
         455             :             !            smaller at the surface.
         456             :             !         2) With the reduction coefficient having a value of 0.2,
         457             :             !            the surface correlation of w & sclr has a value of
         458             :             !            about 0.845.  The correlation is greater if zeta < 0.
         459             :             !            The correlation has a value greater than 1 if
         460             :             !            zeta <= -0.212.
         461             :             !         3) The surface correlations of both rt & sclr and
         462             :             !            thl & sclr are 1.
         463             :             ! Brian Griffin; February 2, 2008.
         464             :             if ( zeta(i) < zero ) then
         465             :
         466             :               sclrprtp(i,1,sclr)  &
         467             :               = reduce_coef  &
         468             :                 * ( wpsclrp_sfc(i,sclr) * wprtp_sfc(i) / ustar(i)**2 ) &
         469             :                 * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
         470             :
         471             :               sclrpthlp(i,1,sclr)  &
         472             :               = reduce_coef  &
         473             :                 * ( wpsclrp_sfc(i,sclr) * wpthlp(i,1) / ustar(i)**2 ) &
         474             :                 * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
         475             :
         476             :               sclrp2(i,1,sclr)  &
         477             :               = reduce_coef   &
         478             :                 * ( wpsclrp_sfc(i,sclr)**2 / ustar(i)**2 ) &
         479             :                 * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
         480             :
         481             :             else
         482             :
         483             :               sclrprtp(i,1,sclr)  &
         484             :               = reduce_coef  &
         485             :                 * four * ( wpsclrp_sfc(i,sclr) * wprtp_sfc(i) / ustar(i)**2 )
         486             :
         487             :               sclrpthlp(i,1,sclr)  &
         488             :               = reduce_coef  &
         489             :                 * four * ( wpsclrp_sfc(i,sclr) * wpthlp(i,1) / ustar(i)**2 )
         490             :
         491             :               sclrp2(i,1,sclr)  &
         492             :               = reduce_coef &
         493             :                 * four * ( wpsclrp_sfc(i,sclr)**2 / ustar(i)**2 )
         494             :
         495             :             endif
         496             :
         497             :           end do
         498             :         enddo ! i = 1, sclr_dim
         499             :         !$acc end parallel loop
         500             :
         501             :       endif
         502             :
         503             :
         504             :     else ! Previous code.
         505             :
         506             :       ! Compute ustar^2
         507             :       !$acc parallel loop gang vector default(present)
         508   149194656 :       do i = 1, ngrdcol
         509   140259600 :         ustar2 = sqrt( upwp_sfc(i) * upwp_sfc(i) + vpwp_sfc(i) * vpwp_sfc(i) )
         510             :
         511             :         ! Compute wstar following Andre et al., 1976
         512   140259600 :         if ( wpthlp(i,1) > zero ) then
         513    93747984 :           if ( .not. l_vary_convect_depth ) then
         514    93747984 :              wstar = ( one/T0 * grav * wpthlp(i,1) * z_const )**(one_third)
         515             :           else
         516           0 :              wstar = ( one/T0 * grav * wpthlp(i,1) * 0.2_core_rknd * depth_pos_wpthlp(i) )**(one_third)
         517             :           end if
         518             :         else
         519             :           wstar = zero
         520             :         endif
         521             :
         522             :         ! Surface friction velocity following Andre et al. 1978
         523   140259600 :         if ( .not. l_vary_convect_depth ) then
         524   140259600 :           uf(i) = sqrt( ustar2 + 0.3_core_rknd * wstar * wstar )
         525             :         else
         526           0 :           uf(i) = sqrt( ustar2 + wstar * wstar )
         527             :         end if
         528             :
         529   149194656 :         uf(i) = max( ufmin, uf(i) )
         530             :       end do
         531             :       !$acc end parallel loop
         532             :
         533             :       ! Compute estimate for surface second order moments
         534             :       !$acc parallel loop gang vector default(present)
         535   149194656 :       do i = 1, ngrdcol
         536   140259600 :         wp2(i,1) = a_const(i) * uf(i)**2
         537   140259600 :         up2(i,1) = up2_sfc_coef(i) * a_const(i) * uf(i)**2  ! From Andre, et al. 1978
         538   149194656 :         vp2(i,1) = up2_sfc_coef(i) * a_const(i) * uf(i)**2  ! "  "
         539             :       end do
         540             :       !$acc end parallel loop
         541             :
         542             :       ! Notes:  1) With "a" having a value of 1.8, the surface correlations of
         543             :       !            both w & rt and w & thl have a value of about 0.878.
         544             :       !         2) The surface correlation of rt & thl is 0.5.
         545             :       ! Brian Griffin; February 2, 2008.
         546             :
         547     8935056 :       if ( .not. l_vary_convect_depth )  then
         548             :         !$acc parallel loop gang vector default(present)
         549   149194656 :         do i = 1, ngrdcol
         550   140259600 :           thlp2(i,1)   = 0.4_core_rknd * a_const(i) * ( wpthlp(i,1) / uf(i) )**2
         551   140259600 :           rtp2(i,1)    = 0.4_core_rknd * a_const(i) * ( wprtp_sfc(i) / uf(i) )**2
         552             :           rtpthlp(i,1) = 0.2_core_rknd * a_const(i) * ( wpthlp(i,1) / uf(i) ) &
         553   149194656 :                                                     * ( wprtp_sfc(i) / uf(i) )
         554             :         end do
         555             :         !$acc end parallel loop
         556             :       else
         557             :         !$acc parallel loop gang vector default(present)
         558           0 :         do i = 1, ngrdcol
         559           0 :           thlp2(i,1)   = ( wpthlp(i,1) / uf(i) )**2 / ( max_mag_correlation_flux**2 * a_const(i) )
         560           0 :           rtp2(i,1)    = ( wprtp_sfc(i) / uf(i) )**2 / ( max_mag_correlation_flux**2 * a_const(i) )
         561           0 :           rtpthlp(i,1) = max_mag_correlation_flux * sqrt( thlp2(i,1) * rtp2(i,1) )
         562             :         end do
         563             :         !$acc end parallel loop
         564             :       end if
         565             :
         566             :       !$acc parallel loop gang vector default(present)
         567   149194656 :       do i = 1, ngrdcol
         568   140259600 :         thlp2(i,1) = max( thl_tol**2, thlp2(i,1) )
         569   149194656 :         rtp2(i,1)  = max( rt_tol**2, rtp2(i,1) )
         570             :       end do
         571             :       !$acc end parallel loop
         572             :
         573             :       ! Add effect of vertical compression of eddies on horizontal gustiness.
         574             :       ! First, ensure that wp2 does not make the correlation
         575             :       !   of (w,thl) or (w,rt) outside (-1,1).
         576             :       ! Perhaps in the future we should also ensure that the correlations
         577             :       !   of (w,u) and (w,v) are not outside (-1,1).
         578             :       !$acc parallel loop gang vector default(present)
         579   149194656 :       do i = 1, ngrdcol
         580   140259600 :         min_wp2_sfc_val(i) &
         581             :             = max( w_tol_sqd, &
         582           0 :                    wprtp_sfc(i)**2 / ( rtp2(i,1) * max_mag_correlation_flux**2 ), &
         583   289454256 :                    wpthlp(i,1)**2 / ( thlp2(i,1) * max_mag_correlation_flux**2 ) )
         584             :       end do
         585             :       !$acc end parallel loop
         586             :
         587             :       !$acc parallel loop gang vector default(present)
         588   149194656 :       do i = 1, ngrdcol
         589   140259600 :         if ( wp2(i,1) - tau_zm(i,1) * lhs_splat_wp2(i,1) * wp2(i,1) &
         590             :              < min_wp2_sfc_val(i) ) then
         591             :         !if ( wp2(i,1) + tau_zm(i,1) * wp2_splat(i,1) < min_wp2_sfc_val(i) ) then
         592             :                            ! splatting correction drives wp2 to overly small values
         593       98548 :           wp2_splat_sfc_correction(i) = -wp2(i,1) + min_wp2_sfc_val(i)
         594       98548 :           wp2(i,1) = min_wp2_sfc_val(i)
         595             :         else
         596   140161052 :          wp2_splat_sfc_correction(i) = tau_zm(i,1) * lhs_splat_wp2(i,1) * wp2(i,1)
         597             :          !wp2_splat_sfc_correction(i) = tau_zm(i,1) * wp2_splat(i,1)
         598   140161052 :          wp2(i,1) = wp2(i,1) + wp2_splat_sfc_correction(i)
         599             :         end if
         600             :
         601   140259600 :         up2(i,1) = up2(i,1) - 0.5_core_rknd * wp2_splat_sfc_correction(i)
         602   149194656 :         vp2(i,1) = vp2(i,1) - 0.5_core_rknd * wp2_splat_sfc_correction(i)
         603             :       end do
         604             :       !$acc end parallel loop
         605             :
         606             :       ! Passive scalars
         607     8935056 :       if ( sclr_dim > 0 ) then
         608             :         !$acc parallel loop gang vector collapse(2) default(present)
         609           0 :         do sclr = 1, sclr_dim
         610           0 :           do i = 1, ngrdcol
         611             :             ! Vince Larson changed coeffs to make correlations between [-1,1].
         612             :             ! 31 Jan 2008
         613             :             ! sclrprtp(i) &
         614             :             ! = a * (wprtp_sfc / uf) * (wpsclrp(i) / uf)
         615             :             ! sclrpthlp(i) &
         616             :             ! = a * (wpthlp / uf) * (wpsclrp(i) / uf)
         617             :             ! sclrp2(i) &
         618             :             ! = sclr_var_coef * a * ( wpsclrp(i) / uf )**2
         619             :             ! Notes:  1) With "a" having a value of 1.8 and "sclr_var_coef"
         620             :             !            having a value of 0.4, the surface correlation of
         621             :             !            w & sclr has a value of about 0.878.
         622             :             !         2) With "sclr_var_coef" having a value of 0.4, the
         623             :             !            surface correlations of both rt & sclr and
         624             :             !            thl & sclr are 0.5.
         625             :             ! Brian Griffin; February 2, 2008.
         626             :
         627             :             ! We use the following if..then's to make sclr_rt and sclr_thl
         628             :             ! close to the actual thlp2/rtp2 at the surface.
         629             :             ! -dschanen 25 Sep 08
         630           0 :             if ( sclr == sclr_idx%iisclr_rt ) then
         631             :               ! If we are trying to emulate rt with the scalar, then we
         632             :               ! use the variance coefficient from above
         633           0 :               sclrprtp(i,1,sclr) = 0.4_core_rknd * a_const(i) &
         634           0 :                                 * ( wprtp_sfc(i) / uf(i) ) * ( wpsclrp_sfc(i,sclr) / uf(i) )
         635             :             else
         636           0 :               sclrprtp(i,1,sclr) = 0.2_core_rknd * a_const(i) &
         637           0 :                                 * ( wprtp_sfc(i) / uf(i) ) * ( wpsclrp_sfc(i,sclr) / uf(i) )
         638             :             endif
         639             :
         640           0 :             if ( sclr == sclr_idx%iisclr_thl ) then
         641             :               ! As above, but for thetal
         642           0 :               sclrpthlp(i,1,sclr) = 0.4_core_rknd * a_const(i) &
         643             :                                  * ( wpthlp(i,1) / uf(i) ) &
         644           0 :                                  * ( wpsclrp_sfc(i,sclr) / uf(i) )
         645             :             else
         646           0 :               sclrpthlp(i,1,sclr) = 0.2_core_rknd * a_const(i) &
         647             :                                  * ( wpthlp(i,1) / uf(i) ) &
         648           0 :                                  * ( wpsclrp_sfc(i,sclr) / uf(i) )
         649             :             endif
         650             :
         651           0 :             sclrp2(i,1,sclr) = sclr_var_coef * a_const(i) &
         652           0 :                                * ( wpsclrp_sfc(i,sclr) / uf(i) )**2
         653             :
         654             :             ! End Vince Larson's change.
         655             :
         656             :           end do
         657             :         enddo ! 1,...sclr_dim
         658             :         !$acc end parallel loop
         659             :
         660             :       endif ! sclr_dim > 0
         661             :
         662             :     endif ! l_andre_1978
         663             :
         664             :     ! Clip wp2 at wp2_max, same as in advance_wp2_wp3
         665             :     !$acc parallel loop gang vector default(present)
         666   149194656 :     do i = 1, ngrdcol
         667   149194656 :       wp2(i,1) = min( wp2(i,1), wp2_max )
         668             :     end do
         669             :     !$acc end parallel loop
         670             :
         671             :     !$acc parallel loop gang vector default(present)
         672   149194656 :     do i = 1, ngrdcol
         673   149194656 :       if ( abs(gr%zm(i,1)-sfc_elevation(i)) > abs(gr%zm(i,1)+sfc_elevation(i))*eps/2 ) then
         674             :
         675             :         ! Variances for cases where the lowest level is not at the surface.
         676             :         ! Eliminate surface effects on lowest level variances.
         677           0 :         wp2(i,1)     = w_tol_sqd
         678           0 :         up2(i,1)     = w_tol_sqd
         679           0 :         vp2(i,1)     = w_tol_sqd
         680           0 :         thlp2(i,1)   = thl_tol**2
         681           0 :         rtp2(i,1)    = rt_tol**2
         682           0 :         rtpthlp(i,1) = 0.0_core_rknd
         683             :
         684             :       end if ! gr%zm(1,1) == sfc_elevation
         685             :     end do
         686             :     !$acc end parallel loop
         687             :
         688     8935056 :     if ( sclr_dim > 0 ) then
         689             :
         690             :       !$acc parallel loop gang vector collapse(2) default(present)
         691           0 :       do sclr = 1, sclr_dim
         692           0 :         do i = 1, ngrdcol
         693           0 :           if ( abs(gr%zm(i,1)-sfc_elevation(i)) > abs(gr%zm(i,1)+sfc_elevation(i))*eps/2 ) then
         694             :
         695             :             ! Variances for cases where the lowest level is not at the surface.
         696             :             ! Eliminate surface effects on lowest level variances.
         697           0 :             sclrp2(i,1,sclr)    = 0.0_core_rknd
         698           0 :             sclrprtp(i,1,sclr)  = 0.0_core_rknd
         699           0 :             sclrpthlp(i,1,sclr) = 0.0_core_rknd
         700             :           end if
         701             :         end do
         702             :       end do
         703             :     end if
         704             :
         705     8935056 :     if ( stats_metadata%l_stats_samp ) then
         706             :
         707             :       !$acc update host( wp2, up2, vp2, thlp2, rtp2, rtpthlp )
         708             :
         709           0 :       do i = 1, ngrdcol
         710             :         call stat_end_update_pt( stats_metadata%ithlp2_sf, 1,    & ! intent(in)
         711           0 :                                  thlp2(i,1) / dt, & ! intent(in)
         712           0 :                                  stats_zm(i) )      ! intent(inout)
         713             :
         714             :         call stat_end_update_pt( stats_metadata%irtp2_sf, 1,     & ! intent(in)
         715             :                                  rtp2(i,1) / dt,  & ! intent(in)
         716           0 :                                  stats_zm(i) )      ! intent(inout)
         717             :
         718             :         call stat_end_update_pt( stats_metadata%irtpthlp_sf, 1,    & ! intent(in)
         719             :                                  rtpthlp(i,1) / dt, & ! intent(in)
         720           0 :                                  stats_zm(i) )        ! intent(inout)
         721             :
         722             :         call stat_end_update_pt( stats_metadata%iup2_sf, 1,    & ! intent(in)
         723             :                                  up2(i,1) / dt, & ! intent(in)
         724           0 :                                  stats_zm(i) )    ! intent(inout)
         725             :
         726             :         call stat_end_update_pt( stats_metadata%ivp2_sf, 1,    & ! intent(in)
         727             :                                  vp2(i,1) / dt, & ! intent(in)
         728           0 :                                  stats_zm(i) )    ! intent(inout)
         729             :
         730             :         call stat_end_update_pt( stats_metadata%iwp2_sf, 1,    & ! intent(in)
         731             :                                  wp2(i,1) / dt, & ! intent(in)
         732           0 :                                  stats_zm(i) )    ! intent(inout)
         733             :       end do
         734             :     end if
         735             :
         736     8935056 :     if ( clubb_at_least_debug_level( 2 ) ) then
         737             :
         738             :       !$acc update host( wp2, up2, vp2, thlp2, rtp2, rtpthlp, &
         739             :       !$acc              upwp_sfc, vpwp_sfc, wpthlp, wprtp_sfc )
         740             :
         741             :       !$acc update host( sclrp2, sclrprtp, sclrpthlp ) if ( sclr_dim > 0 )
         742             :
         743           0 :       do i = 1, ngrdcol
         744           0 :         call sfc_varnce_check( sclr_dim, wp2(i,1), up2(i,1), vp2(i,1),          & ! intent(in)
         745             :                                thlp2(i,1), rtp2(i,1), rtpthlp(i,1),             & ! intent(in)
         746           0 :                                sclrp2(i,1,:), sclrprtp(i,1,:), sclrpthlp(i,1,:) ) ! intent(in)
         747             :       end do
         748             :
         749           0 :       if ( err_code == clubb_fatal_error ) then
         750             :
         751           0 :         write(fstderr,*) "Error in calc_sfc_varnce"
         752           0 :         write(fstderr,*) "Intent(in)"
         753             :
         754           0 :         write(fstderr,*) "upwp_sfc = ", upwp_sfc
         755           0 :         write(fstderr,*) "vpwp_sfc = ", vpwp_sfc
         756           0 :         write(fstderr,*) "wpthlp = ", wpthlp
         757           0 :         write(fstderr,*) "wprtp_sfc = ", wprtp_sfc
         758             :
         759           0 :         if ( sclr_dim > 0 ) then
         760           0 :            write(fstderr,*) "wpsclrp = ", wpsclrp_sfc
         761             :         endif
         762             :
         763           0 :         write(fstderr,*) "Intent(out)"
         764             :
         765           0 :         write(fstderr,*) "wp2 = ", wp2
         766           0 :         write(fstderr,*) "up2 = ", up2
         767           0 :         write(fstderr,*) "vp2 = ", vp2
         768           0 :         write(fstderr,*) "thlp2 = ", thlp2
         769           0 :         write(fstderr,*) "rtp2 = ", rtp2
         770           0 :         write(fstderr,*) "rtpthlp = ", rtpthlp
         771             :
         772           0 :         if ( sclr_dim > 0 ) then
         773           0 :            write(fstderr,*) "sclrp2 = ", sclrp2
         774           0 :            write(fstderr,*) "sclrprtp = ", sclrprtp
         775           0 :            write(fstderr,*) "sclrpthlp = ", sclrpthlp
         776             :         endif
         777             :
         778             :       endif ! err_code == clubb_fatal_error
         779             :
         780             :     endif ! clubb_at_least_debug_level ( 2 )! Update surface stats
         781             :
         782             :     !$acc exit data delete( uf, depth_pos_wpthlp, min_wp2_sfc_val, &
         783             :     !$acc                   um_sfc_sqd, vm_sfc_sqd, usp2_sfc, vsp2_sfc, &
         784             :     !$acc                   ustar, zeta, wp2_splat_sfc_correction )
         785             :
         786     8935056 :     return
         787             :
         788             :   end subroutine calc_sfc_varnce
         789             :
         790             : !===============================================================================
         791             :
         792             : end module sfc_varnce_module

Generated by: LCOV version 1.14