!pbl_dry_contrgradient.f90
!>@brief implements mass-flux contr-gradient parametrization
module pbl_dry_contrgradient
implicit none
type pblContrGradDataType
real, allocatable :: wu(:), &
theta_up(:), &
thetav_up(:), &
ua_up(:), &
va_up(:), &
qv_up(:), &
entr(:)
real :: theta_up0, thetav_up0, qv_up0
real :: ua_up0, va_up0
real :: ua_upH, va_upH
real :: theta_upH, thetav_upH, qv_upH
end type pblContrGradDataType
real, parameter :: c_entr = 0.5 !< dry entrainment parameterization
real, parameter :: b1 = 1.0; !< vertical velocity in updraft Bernulli equation constant
real, parameter :: b2 = 2.0; !< vertical velocity in updraft Bernulli equation constant
public cbl_allocate, cbl_deallocate
public cbl_get_entrainment, cbl_get_up
public cbl_get_wu, cbl_apply_cntrg
contains
subroutine cbl_allocate(cbl, kmax)
integer, intent(in):: kmax
type(pblContrGradDataType), intent(out):: cbl
allocate(cbl%entr(kmax))
allocate(cbl%theta_up(kmax))
allocate(cbl%thetav_up(kmax))
allocate(cbl%qv_up(kmax))
allocate(cbl%ua_up(kmax))
allocate(cbl%va_up(kmax))
allocate(cbl%wu(kmax))
end subroutine cbl_allocate
subroutine cbl_deallocate(cbl)
type(pblContrGradDataType), intent(inout):: cbl
deallocate(cbl%entr)
deallocate(cbl%theta_up)
deallocate(cbl%thetav_up)
deallocate(cbl%qv_up)
deallocate(cbl%ua_up)
deallocate(cbl%va_up)
deallocate(cbl%wu)
end subroutine cbl_deallocate
subroutine cbl_get_entrainment(cbl, bl, fluid, grid)
use scm_state_data, only : stateBLDataType
use pbl_turb_data, only : turbBLDataType
use phys_fluid, only: fluidParamsDataType
use pbl_grid, only : pblgridDataType
implicit none
type(pblContrGradDataType) :: cbl
type(stateBLDataType), intent(inout):: bl
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
integer k
real dz
do k = grid%kmax, 2, -1
dz =(grid%z_cell(k-1) - grid%z_cell(k))
if (grid%z_cell(k) > bl%hpbla_diag) then
cbl%entr(k) = 0.0
else
cbl%entr(k) = c_entr * ( 1/(grid%z_cell(k) + dz) &
+ 1/(bl%hpbla_diag - grid%z_cell(k) + dz));
end if
end do
cbl%entr(1) = 0
end subroutine cbl_get_entrainment
subroutine cbl_get_up(cbl, bl, fluid, grid)
use scm_state_data, only : stateBLDataType
use pbl_turb_data, only : turbBLDataType
use phys_fluid, only: fluidParamsDataType
use pbl_grid, only : pblgridDataType
implicit none
type(pblContrGradDataType) :: cbl
type(stateBLDataType), intent(inout):: bl
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
real phys_beta, wstr, dz
integer k, kmax
kmax = grid%kmax
phys_beta = 0.5 * fluid%g /(bl%theta_v(kmax) + bl%theta_v(kmax-1))
wstr = (phys_beta * (-bl%surf%hs/bl%rho(kmax)) &
* bl%hpbla_diag)**(0.33333)
cbl%thetav_up0 = (-bl%surf%hs/bl%rho(kmax)) / wstr
cbl%theta_up0 = (-bl%surf%hs/bl%rho(kmax)) / wstr
cbl%qv_up0 = (-bl%surf%es/bl%rho(kmax)) / wstr
cbl%ua_up0 = (-bl%surf%cm2u * bl%u(kmax)/bl%rho(kmax)) / wstr
cbl%va_up0 = (-bl%surf%cm2u * bl%v(kmax)/bl%rho(kmax)) / wstr
cbl%thetav_up(kmax) = cbl%thetav_up0
cbl%theta_up(kmax) = cbl%theta_up0
cbl%qv_up(kmax) = cbl%qv_up0
cbl%ua_up(kmax) = cbl%ua_up0
cbl%va_up(kmax) = cbl%va_up0
do k = kmax-1 ,1,-1
if (grid%z_cell(k) > bl%hpbla_diag ) then
cbl%theta_up(k) = 0.0
cbl%thetav_up(k) = 0.0
cbl%qv_up(k) = 0
else
dz = (grid%z_cell(k) - grid%z_cell(k+1))
cbl%theta_up(k) = cbl%theta_up(k + 1) &
- dz * cbl%entr(k+1) * cbl%theta_up(k+1)
cbl%thetav_up(k) = cbl%thetav_up(k + 1) &
- dz * cbl%entr(k+1) * (cbl%thetav_up(k+1))
cbl%qv_up(k) = cbl%qv_up(k + 1) &
- dz * cbl%entr(k+1) * (cbl%qv_up(k+1))
cbl%ua_up(k) = cbl%ua_up(k + 1) &
- dz * cbl%entr(k+1) * (cbl%ua_up(k+1))
cbl%ua_up(k) = cbl%va_up(k + 1) &
- dz * cbl%entr(k+1) * (cbl%va_up(k+1))
end if
end do
do k = kmax ,1,-1
cbl%thetav_up(k) = cbl%thetav_up(k) + bl%theta_v(k)
cbl%theta_up(k) = cbl%theta_up(k) + bl%theta(k)
cbl%qv_up(k) = cbl%qv_up(k) + bl%qv(k)
cbl%ua_up(k) = cbl%ua_up(k) + bl%u(k)
cbl%va_up(k) = cbl%va_up(k) + bl%v(k)
end do
end subroutine cbl_get_up
subroutine cbl_get_wu(cbl, bl, fluid, grid)
use scm_state_data, only : stateBLDataType
use pbl_turb_data, only : turbBLDataType
use phys_fluid, only: fluidParamsDataType
use pbl_grid, only : pblgridDataType
implicit none
type(pblContrGradDataType) :: cbl
type(stateBLDataType), intent(inout):: bl
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
real wu, phys_beta, umod, ustr, rhs, a, b, dz
integer k, kmax
kmax = grid%kmax
cbl%wu(:) = 0.0
umod = sqrt(bl%u(kmax)*bl%u(kmax) + bl%v(kmax)*bl%v(kmax))
ustr = sqrt(bl%surf%cm2u /bl%rho(kmax) * umod)
phys_beta = 0.5 * fluid%g /(bl%theta_v(kmax) + bl%theta_v(kmax-1))
cbl%wu(kmax) = 2.46 * ((ustr)**3 &
- 0.4 * phys_beta * bl%surf%hs &
* (grid%z_cell(k))/bl%rho(kmax))**(1.0/3.0) &
* sqrt(1.0 - grid%z_cell(k) / bl%hpbla_diag)
do k = kmax - 1, 1, -1
dz = (grid%z_cell(K) - grid%z_cell(K+1))
a = 0.0
if (abs(cbl%entr(k)) >= 1.0e-8 .and. abs(cbl%wu(k+1)) > 1.0e-6 ) then
a = b1 * cbl%entr(k) * cbl%wu(k+1) * dz
b = b2 * 9.81 * (dz/cbl%wu(k+1)) &
* (cbl%thetav_up(k) - bl%theta_v(k))/bl%theta_v(k);
rhs = cbl%wu(k+1) - a + b;
if (rhs > 0) then
cbl%wu(k) = rhs;
else
cbl%wu(k) = 0.0;
endif
else
cbl%wu(k) = 0.0;
endif
enddo
end subroutine cbl_get_wu
subroutine cbl_apply_cntrg(cbl, bl, fluid, grid, dt)
use scm_state_data, only : stateBLDataType
use pbl_turb_data, only : turbBLDataType
use phys_fluid, only: fluidParamsDataType
use pbl_grid, only : pblgridDataType
implicit none
type(pblContrGradDataType), intent(inout) :: cbl
type(stateBLDataType), intent(inout):: bl
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
real, intent(in) :: dt
real, dimension(grid%kmax) :: rhs_up, rhs
real a0w0
integer kmax, kpbl, k
kmax = grid%kmax
kpbl = bl%kpbl
a0w0 = 0.15 * cbl%wu(kpbl)
call cbl_get_entrainment(cbl, bl, fluid, grid)
call cbl_get_up(cbl, bl, fluid, grid)
call cbl_get_wu(cbl, bl, fluid, grid)
!apply upwind
rhs_up(:) = 0
do k = kmax-1, 1, -1
if (abs(cbl%wu(k)) > 1.0e-8) then
rhs_up(k) = 0.05 * (cbl%wu(k)) * &
(cbl%theta_up(k+1) - cbl%theta_up(k) + bl%theta(k+1) - bl%theta(k)) &
/ (grid%z_cell(k) - grid%z_cell(k+1))
else
rhs_up(k)=0
end if
end do
bl%theta = bl%theta + dt * rhs_up
rhs_up = 0
do k = kmax-1, 1, -1
if (abs(cbl%wu(k)) > 1.0e-7) then
rhs_up(k) = 0.05 * (cbl%wu(k)) * &
(cbl%qv_up(k+1) - cbl%qv_up(k) &
+ bl%qv(k+1) - bl%qv(k)) / (grid%z_cell(K) -grid%z_cell(K+1))
else
rhs_up(k)=0
end if
end do
bl%qv(:) = bl%qv(:) + dt * rhs_up(:)
end subroutine cbl_apply_cntrg
end module pbl_dry_contrgradient