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

updated pbl

parent beee1a50
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
config-examples/config-gabls1-ex.txt
+ 1
− 1
View file @ c43f1723
......@@ -24,7 +24,7 @@ fluid {
grid {
type = "uniform";
nz = 16;
nz = 64;
h_bot = 0.0;
h_top = 800.0;
}
......
src/diag_pbl.f90
+ 2
− 2
View file @ c43f1723
......@@ -46,7 +46,7 @@ contains
write(*,*) 'here_hpbl'
dz_low = z(kl) - z_surf
hdyn = max(dz_low , 0.5E0 * ustar/cor)
hdyn = max(dz_low , 0.1E0 * ustar/cor)
write(*,*) 'hdyn', hdyn
kpblc = kl
kpbld = kl
......@@ -55,7 +55,7 @@ contains
dz_hdyn = z(k) - z_surf - HDYN
dz_conv = theta_v(k) - theta_v(kl)
if(kpbld.EQ.kl .AND. dz_hdyn.GE.0.E0) kpbld = k
if(kpblc.EQ.kl .AND. dz_conv.GT.0.E0) kpblc = k
if(kpblc.EQ.kl .AND. dz_conv.GT.0.E0) kpblc = k + 1
enddo
kpbl = MIN0(kpblc, kpbld, kl-2)
hpbla = z(kpbl) - z_surf
......
src/pbl_dry_contrgradient.f90
+ 33
− 9
View file @ c43f1723
......@@ -18,7 +18,7 @@ module pbl_dry_contrgradient
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
real, parameter :: b2 = 2.0/3.0; !< vertical velocity in updraft Bernulli equation constant
public cbl_allocate, cbl_deallocate
......@@ -66,7 +66,7 @@ module pbl_dry_contrgradient
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
if (k <= bl%kpbl) then
cbl%entr(k) = 0.0
else
cbl%entr(k) = c_entr * ( 1/(grid%z_cell(k) + dz) &
......@@ -96,8 +96,8 @@ module pbl_dry_contrgradient
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%ua_up0 = 0.0
cbl%va_up0 = 0.0
cbl%thetav_up(kmax) = cbl%thetav_up0
cbl%theta_up(kmax) = cbl%theta_up0
......@@ -106,7 +106,7 @@ module pbl_dry_contrgradient
cbl%va_up(kmax) = cbl%va_up0
do k = kmax-1 ,1,-1
if (grid%z_cell(k) >= bl%hpbla_diag ) then
if (k <= bl%kpbl ) then
cbl%theta_up(k) = 0.0
cbl%thetav_up(k) = 0.0
cbl%qv_up(k) = 0
......@@ -205,8 +205,8 @@ module pbl_dry_contrgradient
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)) &
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
......@@ -218,13 +218,37 @@ module pbl_dry_contrgradient
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))
(cbl%qv_up(+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(:)
!apply upwind
rhs_up(:) = 0
do k = kmax-1, 1, -1
if (abs(cbl%wu(k)) > 1.0e-8) then
rhs_up(k) = 0.00 * (cbl%wu(k)) * &
(cbl%ua_up(k+1) - cbl%ua_up(k) - bl%u(k+1) + bl%u(k)) &
/ (grid%z_cell(k) - grid%z_cell(k+1))
else
rhs_up(k)=0
end if
end do
bl%u(:) = bl%u(:) + dt * rhs_up(:)
rhs_up(:) = 0
do k = kmax-1, 1, -1
if (abs(cbl%wu(k)) > 1.0e-8) then
rhs_up(k) = 0.00 * (cbl%wu(k)) * &
(cbl%va_up(k+1) - cbl%va_up(k) - bl%v(k+1) + bl%v(k)) &
/ (grid%z_cell(k) - grid%z_cell(k+1))
else
rhs_up(k)=0
end if
end do
bl%v(:) = bl%v(:) + dt * rhs_up(:)
end subroutine cbl_apply_cntrg
end module pbl_dry_contrgradient
\ No newline at end of file
src/pbl_fo_turb.f90
+ 42
− 5
View file @ c43f1723
......@@ -8,8 +8,8 @@ module pbl_fo_turb
real, parameter:: yd = 5.e0
real:: ye = yb*yc*sqrt(1.e0/3.e0)
real:: ricr = 2.e0/(3.e0*yd)
real, parameter:: vlinfm = 160.0e0
real:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)
real, parameter:: vlinfm = 60.0e0
real, parameter:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)
real, parameter:: aka = 0.4e0
!real, parameter:: r = 287.05
!real, parameter:: cp = 1.005
......@@ -46,7 +46,8 @@ module pbl_fo_turb
type(turbBLDataType), intent(inout):: turb
if (turb%sf_type == 1) then !> Louis stab functions
call get_fo_coeffs_louis(turb, bl, fluid, grid)
!call get_fo_coeffs_louis(turb, bl, fluid, grid)
call get_fo_coeffs_esau(turb, bl, fluid, grid)
elseif (turb%sf_type == 2) then !> Esau stab functions
! call get_fo_coeffs_esau(turb, bl, fluid, grid)
elseif (turb%sf_type == 3) then !> EFB stab functions
......@@ -91,6 +92,42 @@ module pbl_fo_turb
end do
end subroutine get_fo_coeffs_louis
subroutine get_fo_coeffs_esau(turb, bl, fluid, grid)
use pbl_turb_data, only: turbBLDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
use phys_fluid, only: fluidParamsDataType
use fo_stability_functions
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(inout):: bl
type(pblgridDataType), intent(in)::grid
type(turbBLDataType), intent(inout):: turb
real fnriuv, fnritq,rho
real ocean_flag
integer k, kmax
ocean_flag = real(bl%land_type)
do k = grid%kmax-1, 1, -1
call esau_vit_stab_functions( fnriuv,fnritq,turb%rig(k), turb%l_turb_m(k), &
& grid%z_cell(k), grid%z_cell(k) - grid%z_cell(k+1))
turb%km(k) = fnriuv * turb%l_turb_m(k) * turb%l_turb_m(k) * sqrt(turb%s2(k))
turb%kh(k) = fnritq * turb%l_turb_h(k) * turb%l_turb_h(k) * sqrt(turb%s2(k))
!write(*,*) 'here_coeffs', k, turb%l_turb_m(k), fnriuv, fnritq
end do
! for compliance with old version
bl%vdctq(grid%kmax) = 0.0
bl%vdcuv(grid%kmax) = 0.0
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)
end do
end subroutine get_fo_coeffs_esau
subroutine get_turb_length(turb, bl, fluid, grid, hbl_option)
use pbl_turb_data, only: turbBLDataType
......@@ -139,8 +176,8 @@ module pbl_fo_turb
lh = 1.0
endif
endif
turb%l_turb_m(k) = l_neutral * lm
turb%l_turb_h (k) = l_neutral * lh
turb%l_turb_m(k) = l_neutral * lm / ( 1.0 + l_neutral / vlinfm )
turb%l_turb_h(k) = l_neutral * lh / ( 1.0 + l_neutral / vlinfh )
end do
end subroutine get_turb_length
......
src/pbl_solver.f90
+ 5
− 6
View file @ c43f1723
......@@ -206,7 +206,6 @@ module pbl_solver
if (.not.(allocated(subs_qv))) then
allocate(subs_qv(grid%kmax), source=0.0)
end if
! central differences
do k = kmax-1, 2, -1
rhop = 0.5*(bl%rho(k) + bl%rho(k+1))
......@@ -238,7 +237,7 @@ module pbl_solver
end do
! bottom boundary
write(*,*) 'ws', grid%z_cell(kmax), grid%z_cell(kmax-1)
write(*,*) 'ws', grid%z_cell(kmax), grid%z_cell(kmax-1), w_s(kmax)
k = kmax
subs_u(k) = w_s(k-1) * (bl%u(k)-bl%u(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1))
......@@ -261,10 +260,10 @@ module pbl_solver
!apply
do k=kmax, 1, -1
bl%u(k) = bl%u(k) + subs_u(k) * dt
bl%v(k) = bl%v(k) + subs_v(k) * dt
bl%theta(k) = bl%theta(k) + subs_theta(k) * dt
bl%qv(k) = bl%qv(k) + subs_qv(k) * dt
bl%u(k) = bl%u(k) - subs_u(k) * dt
bl%v(k) = bl%v(k) - subs_v(k) * dt
bl%theta(k) = bl%theta(k) - subs_theta(k) * dt
bl%qv(k) = bl%qv(k) - subs_qv(k) * dt
end do
end subroutine apply_subsidence
end module pbl_solver
\ No newline at end of file
src/tests/bomex.f90 0 → 100644
+ 5
− 0
View file @ c43f1723
! Created by Andrey Debolskiy on 04.04.2025.
program bomex
end program bomex
\ No newline at end of file
src/tests/cbl_exp.f90
+ 3
− 2
View file @ c43f1723
......@@ -299,6 +299,7 @@ program cbl_exp
meteo_cell%dQ = 0.0
meteo_cell%h = grid%z_cell(kmax)
meteo_cell%z0_m = z0
meteo_cell%z0_t = z0/10
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
......@@ -308,14 +309,14 @@ program cbl_exp
tauy = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U * bl_old%v(kmax)
hflx = - bl%rho(kmax) &
* tw_flux0
bl%surf%hs = bl%rho(grid%kmax) * hflx
bl%surf%hs = 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, 2, grid)
call get_coeffs_general(turb, bl, fluid_params, 1, grid)
write(*,*) 'bl%kpbl', bl%kpbl, grid%z_cell(bl%kpbl), sfx_cell%Cm, sfx_cell%Ct, tsurf
call to_file_1d_3var('coeffs1.dat',grid%z_cell, turb%km, turb%km, bl_old%kmax)
!do k=1,grid%kmax
......
src/tests/dycoms2.f90 0 → 100644
+ 5
− 0
View file @ c43f1723
! Created by Andrey Debolskiy on 04.04.2025.
program dycoms2
end program dycoms2
\ No newline at end of file
src/tests/gabls2.f90
+ 11
− 9
View file @ c43f1723
......@@ -248,12 +248,13 @@ program gabls2
tsurf = get_ts_gabls2(time_current)
meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%u(kmax)**2)
meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%v(kmax)**2)
meteo_cell%dT = -tsurf + bl_old%theta(kmax)
meteo_cell%Tsemi = 0.5 * (tsurf + bl_old%theta(kmax))
meteo_cell%dQ = 0.0
meteo_cell%h = grid%z_cell(kmax)
meteo_cell%z0_m = z0
meteo_cell%z0_t= z0 / 10.0
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
......@@ -267,22 +268,22 @@ program gabls2
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
write(*,*) 'surf', bl%surf%hs, turb%ustr, tsurf, bl%theta(kmax)
!write(*,*) 'surf', bl%surf%hs, turb%ustr, tsurf, bl%theta(kmax)
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, grid%z_cell(bl%kpbl), sfx_cell%Cm, sfx_cell%Ct, tsurf
!write(*,*) 'bl%kpbl', bl%kpbl, grid%z_cell(bl%kpbl), sfx_cell%Cm, sfx_cell%Ct, tsurf
call to_file_1d_3var('coeffs1.dat',grid%z_cell, turb%km, turb%km, bl_old%kmax)
!do k=1,grid%kmax
! write(*,*) 'kdiffs', k, bl%vdctq(k), bl%vdcuv(k)
!end do
call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
call cbl_apply_cntrg(cbl, bl, fluid_params, grid, dt)
write(*,*) 'dttheta,', bl%theta- bl_old%theta
!call cbl_apply_cntrg(cbl, bl, fluid_params, grid, dt)
!write(*,*) 'dttheta,', bl%theta - bl_old%theta
call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
write(*,*) 'aftercoriolis,'
!write(*,*) 'aftercoriolis,'
call get_wsub(w_sub, time_current, grid)
!call apply_subsidence(bl, w_sub, fluid_params, grid, dt)
call apply_subsidence(bl, w_sub, fluid_params, grid, dt)
!call to_file_1d_2var('temperature.dat', grid%z_cell, bl_old%U - bl%U, bl_old%kmax)
call scm_data_copy(bl_old,bl)
......@@ -479,10 +480,11 @@ subroutine get_wsub(w, curtime, grid)
if(curtime >= 24.0* 3600.0) then
do k= grid%kmax, 1, -1
if (grid%z_edge(k) <= 1000.0) then
w(k) = -0.005 * (grid%z_edge(k) - 1000.0)
w(k) = - 0.005 * grid%z_edge(k) / 1000.0
else
w(k) = - 0.005
end if
write(*,*) 'ws', k, w(k)
end do
else
w(1:grid%kmax) = 0.0
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment