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