!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