Skip to content
Snippets Groups Projects
Commit 47f3b3a0 authored by Debolskiy Andrey's avatar Debolskiy Andrey :bicyclist_tone5:
Browse files

cntrgradient rework

parent ead9cb39
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 5
− 1
View file @ 47f3b3a0
......@@ -39,6 +39,7 @@ set(lib_files
src/fo_stability_functions.f90
src/pbl_fo_turb.f90
src/pbl_turb_common.f90
src/pbl_dry_contrgradient.f90
src/diag_pbl.f90)
add_library(abl_lib ${lib_files})
......@@ -91,7 +92,10 @@ if (WITH_TESTS)
list(APPEND test_files src/rrtm_interface.f90)
endif()
#gabls1 experiment
add_executable(gabls1 src/tests/gabls1.f90 src/config-utils.f90 src/scm_io_default.f90 src/scm_sfx_data.f90)
add_executable(gabls1 src/tests/gabls1.f90
src/config-utils.f90
src/scm_io_default.f90
src/scm_sfx_data.f90)
target_include_directories(gabls1 PUBLIC ${CMAKE_BINARY_DIR}/modules/)
target_link_libraries( gabls1 PRIVATE abl_lib )
target_link_libraries(gabls1 PRIVATE sfx_lib)
......
config-examples/config-gabls1-ex.txt
+ 1
− 1
View file @ 47f3b3a0
......@@ -8,7 +8,7 @@ fcoriolis = 1.39 * 0.0001;
time{
begin = 0.0;
end = 20.0; #9.0 * 3600.0;
end = 9.0 * 3600.0;
dt = 20.0;
out_dt = 360.0;
}
......
src/pbl_dry_contrgradient.f90 0 → 100644
+ 223
− 0
View file @ 47f3b3a0
!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 allocate_cbl, deallocate_cbl
public get_entrainment, get_up, get_wu, apply_cntrg
contains
subroutine allocate_cbl(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 allocate_cbl
subroutine deallocate_cbl(cbl)
type(pblContrGradDataType), intent(out):: 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 deallocate_cbl
subroutine get_entrainment(cbl, turb, 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(turbBLDataType), intent(in):: turb
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 get_entrainment
subroutine 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 get_up
subroutine 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 get_wu
subroutine 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) :: 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)
!apply upwind
rhs_up(:) = 0
do k = kmax-1, 1, -1
if (abs(cbl%wu(k)) > 1.0e-8) then
rhs_up(k) = 0.1 * (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.1 * (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 apply_cntrg
end module pbl_dry_contrgradient
\ No newline at end of file
src/pbl_fo_turb.f90
+ 0
− 5
View file @ 47f3b3a0
......@@ -74,7 +74,6 @@ module pbl_fo_turb
integer k, kmax
ocean_flag = real(bl%land_type)
write(*,*) 'here'
do k = grid%kmax-1, 1, -1
call louis_stab_functions(fnriuv, fnritq, turb%rig(k), turb%l_turb_m(k), &
& turb%l_turb_m(k), grid%z_cell(k),ocean_flag)
......@@ -84,16 +83,12 @@ module pbl_fo_turb
!write(*,*) 'here_coeffs', k, turb%l_turb_m(k), fnriuv, fnritq
end do
! for compliance with old version
write(*,*) 'here'
bl%vdcuv(grid%kmax) = 0.0
bl%vdcuv(grid%kmax) = 0.0
write(*,*) 'here3.25'
do k = grid%kmax-1,1, -1
bl%vdcuv(k) = 0.5 * (bl%rho(k) + bl%rho(k+1)) * turb%km(k)
bl%vdctq(k) = 0.5 * (bl%rho(k) + bl%rho(k+1)) * turb%kh(k)
write(*,*) bl%vdctq(k), bl%vdcuv(k), turb%km(k)
end do
write(*,*) 'here3.5'
end subroutine get_fo_coeffs_louis
......
src/pbl_solver.f90
+ 18
− 16
View file @ 47f3b3a0
......@@ -44,14 +44,13 @@ module pbl_solver
integer :: k
real :: prgnb(kl)
write(*,*) 'here_diff3.1', bb(ktvd)
! write(*,*) 'here_diff3.1', bb(ktvd)
prgnb(ktvd) = ff(ktvd) / bb(ktvd)
write(*,*) 'here_diff3.1'
!write(*,*) 'here_diff3.1'
do k = ktvd+1, kl
write(*,*) k, ktvd, prgnz(k), ff(k)
!write(*,*) k, ktvd, prgnz(k), ff(k)
prgnb(k) = prgnz(k) * (ff(k) + aa(k) * prgnb(k-1))
end do
write(*,*) 'here_diff3'
y(kl) = prgnb(kl)
do k = kl-1, ktvd, -1
y(k) = prgna(k) * y(k+1) + prgnb(k)
......@@ -76,21 +75,22 @@ module pbl_solver
bb(1:kbltop-1) = 0.0
cc(1:kbltop-1) = 0.0
!top boundary condition: flux = 0
write(*,*) 'here_diff1.25', kbltop
!write(*,*) 'here_diff1.25', kbltop
k = kbltop
dtdz = dt / (grid%dzc(k))
aa(k) = 0
cc(k) = (kdiff(k)/rho(k)) * dtdz / grid%dze(k)
write(*,*) 'here_diff1.5', kdiff(k), kbltop
write(*,*) 'KTVDM', k, aa(k), bb(k), cc(k), rho(k)
bb(k) = 1.0 + cc(k)
!write(*,*) 'here_diff1.5', kdiff(k), kbltop
!write(*,*) 'KTVDM', k, aa(k), bb(k), cc(k), rho(k)
do k = kbltop + 1, grid%kmax -1
dtdz = dt / (grid%dzc(k))
aa(k) = (kdiff(k - 1)/rho(k)) * dtdz / grid%dze(k-1)
cc(k) = (kdiff(k)/rho(k)) * dtdz / grid%dze(k)
bb(k) = 1.0 + aa(k) + cc(k)
write(*,*) 'fill', k, aa(k), bb(k), cc(k), kdiff(k)
!write(*,*) 'fill', k, aa(k), bb(k), cc(k), kdiff(k)
end do
write(*,*) 'here_diff1.75'
!write(*,*) 'here_diff1.75'
!bottom boundary
k = grid%kmax
dtdz = dt / (grid%dzc(k))
......@@ -112,11 +112,11 @@ module pbl_solver
type(pblgridDataType), intent(in) :: grid
real, intent(in):: dt
real, allocatable:: aa(:), bb(:), cc(:), ff(:)
real, save, allocatable:: aa(:), bb(:), cc(:), ff(:)
real, allocatable:: prgna(:), prgnz(:)
integer k, integer, ktop, kmax
kmax = grid%kmax
write(*,*) 'here_diff'
!write(*,*) 'here_diff'
if (.not.(allocated(aa))) then
allocate(aa(grid%kmax), source=0.0)
end if
......@@ -135,7 +135,7 @@ module pbl_solver
if (.not.(allocated(prgnz))) then
allocate(prgnz(grid%kmax), source=0.0)
end if
write(*,*) 'here_diff1', bl%kpbl
!write(*,*) 'here_diff1', bl%kpbl
ktop = bl%kpbl
do k = bl%kpbl-1,1,-1
if(bl%vdcuv(k) > 0.e0) then
......@@ -144,16 +144,16 @@ module pbl_solver
enddo
! fill for temperature and specific humidity
call fill_tridiag(aa, bb, cc, bl%rho, bl%vdctq, ktop, 0.0, grid, dt)
write(*,*) 'here_diff2'
! write(*,*) 'here_diff2'
call factorize_tridiag(ktop, kmax, aa, bb, cc, prgna, prgnz)
do k = ktop,kmax-1
ff(k) = bl%theta(k)
write(*,*) '2', k, aa(k), bb(k), cc(k), ff(k), grid%dze(k)
!write(*,*) '2', k, aa(k), bb(k), cc(k), ff(k), grid%dze(k)
end do
ff(kmax) = bl%theta(kmax) + (dt/kmax) * bl%surf%hs /bl%rho(kmax)
write(*,*) 'here_diff3'
!write(*,*) 'here_diff3'
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%theta)
write(*,*) 'here_diff4'
!write(*,*) 'here_diff4'
do k = ktop,kmax-1
ff(k) = bl%qv(k)
end do
......@@ -174,4 +174,6 @@ module pbl_solver
ff(kmax) = bl%v(kmax)
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%v)
end subroutine solve_diffusion
end module pbl_solver
\ No newline at end of file
src/pbl_turb_common.f90
+ 2
− 4
View file @ 47f3b3a0
......@@ -89,18 +89,16 @@ contains
integer, intent(in):: hbl_option
write(*,*) 'here1'
call get_n2(turb, bl, fluid, grid)
write(*,*) 'here2'
call get_s2(turb, bl, fluid, grid, du2min)
call get_rig(turb, bl, fluid, grid, du2min, ricr)
call get_turb_length(turb, bl, fluid, grid, hbl_option)
if (turb%cl_type == 1) then ! FOM closures
write(*,*) 'here3'
call get_fo_coeffs(turb, bl, fluid, grid)
write(*,*) 'here4'
end if
if (turb%cl_type == 2) then ! KL closures
!call get_kl_coeffs(turb, bl, fluid, grid)
......
src/tests/gabls1.f90
+ 10
− 7
View file @ 47f3b3a0
......@@ -258,8 +258,9 @@ program gabls1
!finish updating grid
call get_sigma_from_z_theta( grid, bl, fluid_params%p0 * 100.0, fluid_params)
call theta2ta(bl_old%temp, bl_old%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call theta2ta(bl%temp, bl%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
call scm_data_copy(bl,bl_old)
call scm_data_copy(bl_old,bl)
!Initialize turbulent closure
call turb_data_allocate(turb, kmax)
......@@ -298,14 +299,15 @@ program gabls1
bl%surf%hs = bl%rho(grid%kmax) * hflx
bl%surf%es = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U *meteo_cell%dQ
turb%ustr = sfx_cell%Cm * meteo_cell%U
umst = turb%ustr
call get_density_theta_vector( bl, fluid_params, grid)
call get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
call get_coeffs_general(turb, bl, fluid_params, 1, grid)
write(*,*) 'bl%kpbl', bl%kpbl
call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
!call add_coriolis(bl, ug, vg, f_cor, dt, grid)
call scm_data_copy(bl,bl_old)
call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
call scm_data_copy(bl_old,bl)
time_current = time_current + dt
......@@ -389,18 +391,19 @@ end subroutine init_state
subroutine add_coriolis(bl, ug, vg, f, dt, grid)
subroutine add_coriolis(bl, bl_old, ug, vg, f, dt, grid)
use scm_state_data
use pbl_grid, only: pblgridDataType
implicit none
real, intent(in) :: ug, vg, f, dt
type(stateBLDataType), intent(out):: bl
type(stateBLDataType), intent(inout):: bl
type(stateBLDataType), intent(in):: bl_old
type(pblgridDataType), intent(in):: grid
integer k
do k = 1, grid%kmax
bl%v(k) = bl%v(k) - dt * f * (bl%u(k) - ug)
bl%u(k) = bl%u(k) + dt * f * ( bl%v(k) - vg)
bl%v(k) = bl_old%v(k) - dt * f * (bl_old%u(k) - ug)
bl%u(k) = bl_old%u(k) + dt * f * ( bl_old%v(k) - vg)
end do
end subroutine add_coriolis
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment