LCOV - code coverage report
Current view: | top level - physics/cam - pbl_utils.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 35 | 92 | 38.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 6 | 10 | 60.0 %

              Line data    Source code
           1             : module pbl_utils
           2             : !-----------------------------------------------------------------------!
           3             : ! Module to hold PBL-related subprograms that may be used with multiple !
           4             : ! different vertical diffusion schemes.                                 !
           5             : !                                                                       !
           6             : ! Public subroutines:                                                   !
           7             : !
           8             : !     calc_obklen                                                       !
           9             : !                                                                       !
          10             : !------------------ History --------------------------------------------!
          11             : ! Created: Apr. 2012, by S. Santos                                      !
          12             : !-----------------------------------------------------------------------!
          13             :
          14             : use shr_kind_mod, only: r8 => shr_kind_r8
          15             :
          16             : implicit none
          17             : private
          18             :
          19             : ! Public Procedures
          20             : !----------------------------------------------------------------------!
          21             : ! Excepting the initialization procedure, these are elemental
          22             : ! procedures, so they can accept scalars or any dimension of array as
          23             : ! arguments, as long as all arguments have the same number of
          24             : ! elements.
          25             : public pbl_utils_init
          26             : public calc_ustar
          27             : public calc_obklen
          28             : public virtem
          29             : public compute_radf
          30             : public austausch_atm, austausch_atm_free
          31             :
          32             : real(r8), parameter :: ustar_min = 0.01_r8
          33             :
          34             : real(r8) :: g         ! acceleration of gravity
          35             : real(r8) :: vk        ! Von Karman's constant
          36             : real(r8) :: cpair     ! specific heat of dry air
          37             : real(r8) :: rair      ! gas constant for dry air
          38             : real(r8) :: zvir      ! rh2o/rair - 1
          39             :
          40             :
          41             : !------------------------------------------------------------------------!
          42             : ! Purpose: Compilers aren't creating optimized vector versions of        !
          43             : !          elemental routines, so we'll explicitly create them and bind  !
          44             : !          them via an interface for transparent use                     !
          45             : !------------------------------------------------------------------------!
          46             : interface calc_ustar
          47             :   module procedure calc_ustar_scalar
          48             :   module procedure calc_ustar_vector
          49             : end interface
          50             :
          51             : interface calc_obklen
          52             :   module procedure calc_obklen_scalar
          53             :   module procedure calc_obklen_vector
          54             : end interface
          55             :
          56             : interface virtem
          57             :   module procedure virtem_vector1D
          58             :   module procedure virtem_vector2D  ! Used in hb_diff.F90
          59             : end interface
          60             :
          61             :
          62             :
          63             : contains
          64             :
          65        1536 : subroutine pbl_utils_init(g_in,vk_in,cpair_in,rair_in,zvir_in)
          66             :
          67             :   !-----------------------------------------------------------------------!
          68             :   ! Purpose: Set constants to be used in calls to later functions         !
          69             :   !-----------------------------------------------------------------------!
          70             :
          71             :   real(r8), intent(in) :: g_in       ! acceleration of gravity
          72             :   real(r8), intent(in) :: vk_in      ! Von Karman's constant
          73             :   real(r8), intent(in) :: cpair_in   ! specific heat of dry air
          74             :   real(r8), intent(in) :: rair_in    ! gas constant for dry air
          75             :   real(r8), intent(in) :: zvir_in    ! rh2o/rair - 1
          76             :
          77        1536 :   g = g_in
          78        1536 :   vk = vk_in
          79        1536 :   cpair = cpair_in
          80        1536 :   rair = rair_in
          81        1536 :   zvir = zvir_in
          82             :
          83        1536 : end subroutine pbl_utils_init
          84             :
          85    70129800 : subroutine calc_ustar_scalar( t,    pmid, taux, tauy, &
          86             :                                  rrho, ustar)
          87             :
          88             :   !-----------------------------------------------------------------------!
          89             :   ! Purpose: Calculate ustar and bottom level density (necessary for      !
          90             :   !  Obukhov length calculation).                                         !
          91             :   !-----------------------------------------------------------------------!
          92             :
          93             :   real(r8), intent(in) :: t         ! surface temperature
          94             :   real(r8), intent(in) :: pmid      ! midpoint pressure (bottom level)
          95             :   real(r8), intent(in) :: taux      ! surface u stress [N/m2]
          96             :   real(r8), intent(in) :: tauy      ! surface v stress [N/m2]
          97             :
          98             :   real(r8), intent(out) :: rrho     ! 1./bottom level density
          99             :   real(r8), intent(out) :: ustar    ! surface friction velocity [m/s]
         100             :
         101    70129800 :   rrho = rair * t / pmid
         102    70129800 :   ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min )
         103             :
         104    70129800 : end subroutine calc_ustar_scalar
         105             :
         106     5956704 : subroutine calc_ustar_vector(n, t, pmid, taux, tauy, &
         107     5956704 :                                  rrho, ustar)
         108             :
         109             :   !-----------------------------------------------------------------------!
         110             :   ! Purpose: Calculate ustar and bottom level density (necessary for      !
         111             :   !  Obukhov length calculation).                                         !
         112             :   !-----------------------------------------------------------------------!
         113             :   integer, intent(in) :: n             ! Length of vectors
         114             :
         115             :   real(r8), intent(in) :: t(n)         ! surface temperature
         116             :   real(r8), intent(in) :: pmid(n)      ! midpoint pressure (bottom level)
         117             :   real(r8), intent(in) :: taux(n)      ! surface u stress [N/m2]
         118             :   real(r8), intent(in) :: tauy(n)      ! surface v stress [N/m2]
         119             :
         120             :
         121             :   real(r8), intent(out) :: rrho(n)     ! 1./bottom level density
         122             :   real(r8), intent(out) :: ustar(n)    ! surface friction velocity [m/s]
         123             :
         124             :
         125    99463104 :   rrho = rair * t / pmid
         126    99463104 :   ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min )
         127             :
         128     5956704 : end subroutine calc_ustar_vector
         129             :
         130           0 : subroutine calc_obklen_scalar( ths,  thvs, qflx, shflx, rrho, ustar, &
         131             :                                   khfs, kqfs, kbfs, obklen)
         132             :
         133             :   !-----------------------------------------------------------------------!
         134             :   ! Purpose: Calculate Obukhov length and kinematic fluxes.               !
         135             :   !-----------------------------------------------------------------------!
         136             :
         137             :   real(r8), intent(in)  :: ths           ! potential temperature at surface [K]
         138             :   real(r8), intent(in)  :: thvs          ! virtual potential temperature at surface
         139             :   real(r8), intent(in)  :: qflx          ! water vapor flux (kg/m2/s)
         140             :   real(r8), intent(in)  :: shflx         ! surface heat flux (W/m2)
         141             :
         142             :   real(r8), intent(in)  :: rrho          ! 1./bottom level density [ m3/kg ]
         143             :   real(r8), intent(in)  :: ustar         ! Surface friction velocity [ m/s ]
         144             :
         145             :   real(r8), intent(out) :: khfs          ! sfc kinematic heat flux [mK/s]
         146             :   real(r8), intent(out) :: kqfs          ! sfc kinematic water vapor flux [m/s]
         147             :   real(r8), intent(out) :: kbfs          ! sfc kinematic buoyancy flux [m^2/s^3]
         148             :   real(r8), intent(out) :: obklen        ! Obukhov length
         149             :
         150             :   ! Need kinematic fluxes for Obukhov:
         151           0 :   khfs = shflx*rrho/cpair
         152           0 :   kqfs = qflx*rrho
         153           0 :   kbfs = khfs + zvir*ths*kqfs
         154             :
         155             :   ! Compute Obukhov length:
         156           0 :   obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs)))
         157             :
         158           0 : end subroutine calc_obklen_scalar
         159             :
         160     5956704 : subroutine calc_obklen_vector(n, ths,  thvs, qflx, shflx, rrho, ustar, &
         161     5956704 :                                   khfs, kqfs, kbfs, obklen)
         162             :
         163             :   !-----------------------------------------------------------------------!
         164             :   ! Purpose: Calculate Obukhov length and kinematic fluxes.               !
         165             :   !-----------------------------------------------------------------------!
         166             :   integer, intent(in) :: n                  ! Length of vectors
         167             :
         168             :   real(r8), intent(in)  :: ths(n)           ! potential temperature at surface [K]
         169             :   real(r8), intent(in)  :: thvs(n)          ! virtual potential temperature at surface
         170             :   real(r8), intent(in)  :: qflx(n)          ! water vapor flux (kg/m2/s)
         171             :   real(r8), intent(in)  :: shflx(n)         ! surface heat flux (W/m2)
         172             :
         173             :   real(r8), intent(in)  :: rrho(n)          ! 1./bottom level density [ m3/kg ]
         174             :   real(r8), intent(in)  :: ustar(n)         ! Surface friction velocity [ m/s ]
         175             :
         176             :   real(r8), intent(out) :: khfs(n)          ! sfc kinematic heat flux [mK/s]
         177             :   real(r8), intent(out) :: kqfs(n)          ! sfc kinematic water vapor flux [m/s]
         178             :   real(r8), intent(out) :: kbfs(n)          ! sfc kinematic buoyancy flux [m^2/s^3]
         179             :   real(r8), intent(out) :: obklen(n)        ! Obukhov length
         180             :
         181             :
         182             :   ! Need kinematic fluxes for Obukhov:
         183    99463104 :   khfs = shflx*rrho/cpair
         184    99463104 :   kqfs = qflx*rrho
         185    99463104 :   kbfs = khfs + zvir*ths*kqfs
         186             :
         187             :   ! Compute Obukhov length:
         188    99463104 :   obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs)))
         189             :
         190     5956704 : end subroutine calc_obklen_vector
         191             :
         192           0 : subroutine virtem_vector1D(n, t,q, virtem)
         193             :
         194             :   !-----------------------------------------------------------------------!
         195             :   ! Purpose: Calculate virtual temperature from temperature and specific  !
         196             :   !  humidity.                                                            !
         197             :   !-----------------------------------------------------------------------!
         198             :
         199             :   integer,  intent(in) :: n              ! vector length
         200             :
         201             :   real(r8), intent(in) :: t(n), q(n)
         202             :   real(r8), intent(out):: virtem(n)
         203             :
         204           0 :   virtem = t * (1.0_r8 + zvir*q)
         205             :
         206           0 : end subroutine virtem_vector1D
         207             :
         208     1489176 : subroutine virtem_vector2D(n, m, t, q, virtem)
         209             :
         210             :   !-----------------------------------------------------------------------!
         211             :   ! Purpose: Calculate virtual temperature from temperature and specific  !
         212             :   !  humidity.                                                            !
         213             :   !-----------------------------------------------------------------------!
         214             :
         215             :   integer,  intent(in) :: n, m            ! vector lengths
         216             :
         217             :   real(r8), intent(in) :: t(n,m), q(n,m)
         218             :   real(r8), intent(out):: virtem(n,m)
         219             :
         220  2314006344 :   virtem = t * (1.0_r8 + zvir*q)
         221             :
         222     1489176 : end subroutine virtem_vector2D
         223             :
         224             :
         225           0 : subroutine compute_radf( choice_radf, i, pcols, pver, ncvmax, ncvfin, ktop, qmin, &
         226           0 :                          ql, pi, qrlw, g, cldeff, zi, chs, lwp_CL, opt_depth_CL,  &
         227           0 :                          radinvfrac_CL, radf_CL )
         228             :   ! -------------------------------------------------------------------------- !
         229             :   ! Purpose:                                                                   !
         230             :   ! Calculate cloud-top radiative cooling contribution to buoyancy production. !
         231             :   ! Here,  'radf' [m2/s3] is additional buoyancy flux at the CL top interface  !
         232             :   ! associated with cloud-top LW cooling being mainly concentrated near the CL !
         233             :   ! top interface ( just below CL top interface ).  Contribution of SW heating !
         234             :   ! within the cloud is not included in this radiative buoyancy production     !
         235             :   ! since SW heating is more broadly distributed throughout the CL top layer.  !
         236             :   ! -------------------------------------------------------------------------- !
         237             :
         238             :   !-----------------!
         239             :   ! Input variables !
         240             :   !-----------------!
         241             :   character(len=6), intent(in) :: choice_radf  ! Method for calculating radf
         242             :   integer,  intent(in)  :: i                   ! Index of current column
         243             :   integer,  intent(in)  :: pcols               ! Number of atmospheric columns
         244             :   integer,  intent(in)  :: pver                ! Number of atmospheric layers
         245             :   integer,  intent(in)  :: ncvmax              ! Max numbers of CLs (perhaps equal to pver)
         246             :   integer,  intent(in)  :: ncvfin(pcols)       ! Total number of CL in column
         247             :   integer,  intent(in)  :: ktop(pcols, ncvmax) ! ktop for current column
         248             :   real(r8), intent(in)  :: qmin                ! Minimum grid-mean LWC counted as clouds [kg/kg]
         249             :   real(r8), intent(in)  :: ql(pcols, pver)     ! Liquid water specific humidity [ kg/kg ]
         250             :   real(r8), intent(in)  :: pi(pcols, pver+1)   ! Interface pressures [ Pa ]
         251             :   real(r8), intent(in)  :: qrlw(pcols, pver)   ! Input grid-mean LW heating rate : [ K/s ] * cpair * dp = [ W/kg*Pa ]
         252             :   real(r8), intent(in)  :: g                   ! Gravitational acceleration
         253             :   real(r8), intent(in)  :: cldeff(pcols,pver)  ! Effective Cloud Fraction [fraction]
         254             :   real(r8), intent(in)  :: zi(pcols, pver+1)   ! Interface heights [ m ]
         255             :   real(r8), intent(in)  :: chs(pcols, pver+1)  ! Buoyancy coeffi. saturated sl (heat) coef. at all interfaces.
         256             :
         257             :   !------------------!
         258             :   ! Output variables !
         259             :   !------------------!
         260             :   real(r8), intent(out) :: lwp_CL(ncvmax)         ! LWP in the CL top layer [ kg/m2 ]
         261             :   real(r8), intent(out) :: opt_depth_CL(ncvmax)   ! Optical depth of the CL top layer
         262             :   real(r8), intent(out) :: radinvfrac_CL(ncvmax)  ! Fraction of LW radiative cooling confined in the top portion of CL
         263             :   real(r8), intent(out) :: radf_CL(ncvmax)        ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
         264             :
         265             :   !-----------------!
         266             :   ! Local variables !
         267             :   !-----------------!
         268             :   integer :: kt, ncv
         269             :   real(r8) :: lwp, opt_depth, radinvfrac, radf
         270             :
         271             :
         272             :   !-----------------!
         273             :   ! Begin main code !
         274             :   !-----------------!
         275           0 :   lwp_CL        = 0._r8
         276           0 :   opt_depth_CL  = 0._r8
         277           0 :   radinvfrac_CL = 0._r8
         278           0 :   radf_CL       = 0._r8
         279             :
         280             :   ! ---------------------------------------- !
         281             :   ! Perform do loop for individual CL regime !
         282             :   ! ---------------------------------------- !
         283           0 :   do ncv = 1, ncvfin(i)
         284           0 :     kt = ktop(i,ncv)
         285             :     !-----------------------------------------------------!
         286             :     ! Compute radf for each CL regime and for each column !
         287             :     !-----------------------------------------------------!
         288           0 :     if( choice_radf .eq. 'orig' ) then
         289           0 :       if( ql(i,kt) .gt. qmin .and. ql(i,kt-1) .lt. qmin ) then
         290           0 :         lwp       = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
         291           0 :         opt_depth = 156._r8 * lwp  ! Estimated LW optical depth in the CL top layer
         292             :         ! Approximate LW cooling fraction concentrated at the inversion by using
         293             :         ! polynomial approx to exact formula 1-2/opt_depth+2/(exp(opt_depth)-1))
         294             :
         295           0 :         radinvfrac  = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
         296           0 :         radf        = qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) ! Cp*radiative cooling = [ W/kg ]
         297           0 :         radf        = max( radinvfrac * radf * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 ) * chs(i,kt)
         298             :         ! We can disable cloud LW cooling contribution to turbulence by uncommenting:
         299             :         ! radf = 0._r8
         300             :       end if
         301             :
         302           0 :     elseif( choice_radf .eq. 'ramp' ) then
         303             :
         304           0 :       lwp         = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
         305           0 :       opt_depth   = 156._r8 * lwp  ! Estimated LW optical depth in the CL top layer
         306           0 :       radinvfrac  = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
         307           0 :       radinvfrac  = max(cldeff(i,kt)-cldeff(i,kt-1),0._r8) * radinvfrac
         308           0 :       radf        = qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) ! Cp*radiative cooling [W/kg]
         309           0 :       radf        = max( radinvfrac * radf * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 ) * chs(i,kt)
         310             :
         311           0 :     elseif( choice_radf .eq. 'maxi' ) then
         312             :
         313             :       ! Radiative flux divergence both in 'kt' and 'kt-1' layers are included
         314             :       ! 1. From 'kt' layer
         315           0 :         lwp         = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
         316           0 :         opt_depth   = 156._r8 * lwp  ! Estimated LW optical depth in the CL top layer
         317           0 :         radinvfrac  = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
         318           0 :         radf        = max( radinvfrac * qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 )
         319             :       ! 2. From 'kt-1' layer and add the contribution from 'kt' layer
         320           0 :         lwp         = ql(i,kt-1) * ( pi(i,kt) - pi(i,kt-1) ) / g
         321           0 :         opt_depth   = 156._r8 * lwp  ! Estimated LW optical depth in the CL top layer
         322           0 :         radinvfrac  = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth) + opt_depth**2 )
         323           0 :         radf        = radf + max( radinvfrac * qrlw(i,kt-1) / ( pi(i,kt-1) - pi(i,kt) ) * ( zi(i,kt-1) - zi(i,kt) ), 0._r8 )
         324           0 :         radf        = max( radf, 0._r8 ) * chs(i,kt)
         325             :
         326             :     endif
         327             :
         328           0 :     lwp_CL(ncv)        = lwp
         329           0 :     opt_depth_CL(ncv)  = opt_depth
         330           0 :     radinvfrac_CL(ncv) = radinvfrac
         331           0 :     radf_CL(ncv)       = radf
         332             :   end do ! ncv = 1, ncvfin(i)
         333           0 : end subroutine compute_radf
         334             :
         335           0 : subroutine austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
         336             :
         337             :   !---------------------------------------------------------------------- !
         338             :   !                                                                       !
         339             :   ! Purpose: Computes exchange coefficients for free turbulent flows.     !
         340             :   !                                                                       !
         341             :   ! Method:                                                               !
         342             :   !                                                                       !
         343             :   ! The free atmosphere diffusivities are based on standard mixing length !
         344             :   ! forms for the neutral diffusivity multiplied by functns of Richardson !
         345             :   ! number. K = l^2 * |dV/dz| * f(Ri). The same functions are used for    !
         346             :   ! momentum, potential temperature, and constitutents.                   !
         347             :   !                                                                       !
         348             :   ! The stable Richardson num function (Ri>0) is taken from Holtslag and  !
         349             :   ! Beljaars (1989), ECMWF proceedings. f = 1 / (1 + 10*Ri*(1 + 8*Ri))    !
         350             :   ! The unstable Richardson number function (Ri<0) is taken from  CCM1.   !
         351             :   ! f = sqrt(1 - 18*Ri)                                                   !
         352             :   !                                                                       !
         353             :   ! Author: B. Stevens (rewrite, August 2000)                             !
         354             :   !                                                                       !
         355             :   !---------------------------------------------------------------------- !
         356             :
         357             :   ! --------------- !
         358             :   ! Input arguments !
         359             :   ! --------------- !
         360             :
         361             :   integer,  intent(in)  :: pcols                ! Atmospheric columns dimension size
         362             :   integer,  intent(in)  :: ncol                 ! Number of atmospheric columns
         363             :   integer,  intent(in)  :: pver                 ! Number of atmospheric layers
         364             :   integer,  intent(in)  :: ntop                 ! Top layer for calculation
         365             :   integer,  intent(in)  :: nbot                 ! Bottom layer for calculation
         366             :
         367             :   real(r8), intent(in)  :: ml2(pver+1)          ! Mixing lengths squared
         368             :   real(r8), intent(in)  :: s2(pcols,pver)       ! Shear squared
         369             :   real(r8), intent(in)  :: ri(pcols,pver)       ! Richardson no
         370             :
         371             :   ! ---------------- !
         372             :   ! Output arguments !
         373             :   ! ---------------- !
         374             :
         375             :   real(r8), intent(out) :: kvf(pcols,pver+1)    ! Eddy diffusivity for heat and tracers
         376             :
         377             :   ! --------------- !
         378             :   ! Local Variables !
         379             :   ! --------------- !
         380             :
         381             :   real(r8)              :: fofri                ! f(ri)
         382             :   real(r8)              :: kvn                  ! Neutral Kv
         383             :
         384             :   integer               :: i                    ! Longitude index
         385             :   integer               :: k                    ! Vertical index
         386             :
         387             :   real(r8), parameter :: zkmin =  0.01_r8       ! Minimum kneutral*f(ri).
         388             :
         389             :   ! ----------------------- !
         390             :   ! Main Computation Begins !
         391             :   ! ----------------------- !
         392             :
         393           0 :   kvf(:ncol,:)           = 0.0_r8
         394             :
         395             :   ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
         396             :
         397           0 :   do k = ntop, nbot - 1
         398           0 :      do i = 1, ncol
         399           0 :         if( ri(i,k) < 0.0_r8 ) then
         400           0 :            fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
         401             :         else
         402           0 :            fofri = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * ri(i,k) * ( 1.0_r8 + 8.0_r8 * ri(i,k) ) )
         403             :         end if
         404           0 :         kvn = ml2(k) * sqrt(s2(i,k))
         405           0 :         kvf(i,k+1) = max( zkmin, kvn * fofri )
         406             :      end do
         407             :   end do
         408             :
         409           0 : end subroutine austausch_atm
         410             :
         411     1489176 : subroutine austausch_atm_free(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
         412             :
         413             :   !---------------------------------------------------------------------- !
         414             :   !                                                                       !
         415             :   ! same as austausch_atm but only mixing for Ri<0                        !
         416             :   ! i.e. no background mixing and mixing for Ri>0                         !
         417             :   !                                                                       !
         418             :   !---------------------------------------------------------------------- !
         419             :
         420             :   ! --------------- !
         421             :   ! Input arguments !
         422             :   ! --------------- !
         423             :
         424             :   integer,  intent(in)  :: pcols                ! Atmospheric columns dimension size
         425             :   integer,  intent(in)  :: ncol                 ! Number of atmospheric columns
         426             :   integer,  intent(in)  :: pver                 ! Number of atmospheric layers
         427             :   integer,  intent(in)  :: ntop                 ! Top layer for calculation
         428             :   integer,  intent(in)  :: nbot                 ! Bottom layer for calculation
         429             :
         430             :   real(r8), intent(in)  :: ml2(pver+1)          ! Mixing lengths squared
         431             :   real(r8), intent(in)  :: s2(pcols,pver)       ! Shear squared
         432             :   real(r8), intent(in)  :: ri(pcols,pver)       ! Richardson no
         433             :
         434             :   ! ---------------- !
         435             :   ! Output arguments !
         436             :   ! ---------------- !
         437             :
         438             :   real(r8), intent(out) :: kvf(pcols,pver+1)    ! Eddy diffusivity for heat and tracers
         439             :
         440             :   ! --------------- !
         441             :   ! Local Variables !
         442             :   ! --------------- !
         443             :
         444             :   real(r8)              :: fofri                ! f(ri)
         445             :   real(r8)              :: kvn                  ! Neutral Kv
         446             :
         447             :   integer               :: i                    ! Longitude index
         448             :   integer               :: k                    ! Vertical index
         449             :
         450             :   ! ----------------------- !
         451             :   ! Main Computation Begins !
         452             :   ! ----------------------- !
         453             :
         454  2338872120 :   kvf(:ncol,:)           = 0.0_r8
         455             :   ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
         456   138493368 :   do k = ntop, nbot - 1
         457  2289140568 :      do i = 1, ncol
         458  2150647200 :         if( ri(i,k) < 0.0_r8 ) then
         459    47375959 :            fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
         460             :         else
         461             :            fofri = 0.0_r8
         462             :         end if
         463  2150647200 :         kvn = ml2(k) * sqrt(s2(i,k))
         464  2287651392 :         kvf(i,k+1) = kvn * fofri
         465             :      end do
         466             :   end do
         467     1489176 : end subroutine austausch_atm_free
         468             :
         469             : end module pbl_utils

Generated by: LCOV version 1.14