From c43f17234e0d38505186cfbf2ef1eadea87b1807 Mon Sep 17 00:00:00 2001
From: Andrey Debolskiy <and.debol@gmail.com>
Date: Fri, 11 Apr 2025 10:42:44 +0300
Subject: [PATCH] updated pbl
---
config-examples/config-gabls1-ex.txt | 2 +-
src/diag_pbl.f90 | 4 +--
src/pbl_dry_contrgradient.f90 | 42 +++++++++++++++++++------
src/pbl_fo_turb.f90 | 47 +++++++++++++++++++++++++---
src/pbl_solver.f90 | 11 +++----
src/tests/bomex.f90 | 5 +++
src/tests/cbl_exp.f90 | 5 +--
src/tests/dycoms2.f90 | 5 +++
src/tests/gabls2.f90 | 20 ++++++------
9 files changed, 107 insertions(+), 34 deletions(-)
create mode 100644 src/tests/bomex.f90
create mode 100644 src/tests/dycoms2.f90
diff --git a/config-examples/config-gabls1-ex.txt b/config-examples/config-gabls1-ex.txt
index 3c6a9bd..b677501 100644
--- a/config-examples/config-gabls1-ex.txt
+++ b/config-examples/config-gabls1-ex.txt
@@ -24,7 +24,7 @@ fluid {
grid {
type = "uniform";
- nz = 16;
+ nz = 64;
h_bot = 0.0;
h_top = 800.0;
}
diff --git a/src/diag_pbl.f90 b/src/diag_pbl.f90
index dbabc0e..2bef69b 100644
--- a/src/diag_pbl.f90
+++ b/src/diag_pbl.f90
@@ -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
diff --git a/src/pbl_dry_contrgradient.f90 b/src/pbl_dry_contrgradient.f90
index 42bc01f..4afbbcc 100644
--- a/src/pbl_dry_contrgradient.f90
+++ b/src/pbl_dry_contrgradient.f90
@@ -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
diff --git a/src/pbl_fo_turb.f90 b/src/pbl_fo_turb.f90
index ab33aff..0ed5830 100644
--- a/src/pbl_fo_turb.f90
+++ b/src/pbl_fo_turb.f90
@@ -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
diff --git a/src/pbl_solver.f90 b/src/pbl_solver.f90
index f7b99a2..ae92ec2 100644
--- a/src/pbl_solver.f90
+++ b/src/pbl_solver.f90
@@ -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
diff --git a/src/tests/bomex.f90 b/src/tests/bomex.f90
new file mode 100644
index 0000000..b06eabe
--- /dev/null
+++ b/src/tests/bomex.f90
@@ -0,0 +1,5 @@
+! Created by Andrey Debolskiy on 04.04.2025.
+
+program bomex
+
+end program bomex
\ No newline at end of file
diff --git a/src/tests/cbl_exp.f90 b/src/tests/cbl_exp.f90
index 9ea944c..824cc96 100644
--- a/src/tests/cbl_exp.f90
+++ b/src/tests/cbl_exp.f90
@@ -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
diff --git a/src/tests/dycoms2.f90 b/src/tests/dycoms2.f90
new file mode 100644
index 0000000..6e0ed5b
--- /dev/null
+++ b/src/tests/dycoms2.f90
@@ -0,0 +1,5 @@
+! Created by Andrey Debolskiy on 04.04.2025.
+
+program dycoms2
+
+end program dycoms2
\ No newline at end of file
diff --git a/src/tests/gabls2.f90 b/src/tests/gabls2.f90
index 90739d3..4ad2edf 100644
--- a/src/tests/gabls2.f90
+++ b/src/tests/gabls2.f90
@@ -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
+ w(k) = - 0.005
end if
+ write(*,*) 'ws', k, w(k)
end do
else
w(1:grid%kmax) = 0.0
--
GitLab