diff --git a/CMakeLists.txt b/CMakeLists.txt
index 46f863338fa9d31da264fa93b2f1089096fc478a..61a7611509808b9558e2afe9f6022a7049906967 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -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)
diff --git a/config-examples/config-gabls1-ex.txt b/config-examples/config-gabls1-ex.txt
index 3b20be4fe4cbfc8b49977d133a5674d2faa92c01..ed577962618120bdb6527972b1ee955efca4b210 100644
--- a/config-examples/config-gabls1-ex.txt
+++ b/config-examples/config-gabls1-ex.txt
@@ -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;
}
diff --git a/src/pbl_dry_contrgradient.f90 b/src/pbl_dry_contrgradient.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ff89bbee78b68fc368d2415485f2738733e6064d
--- /dev/null
+++ b/src/pbl_dry_contrgradient.f90
@@ -0,0 +1,223 @@
+!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
diff --git a/src/pbl_fo_turb.f90 b/src/pbl_fo_turb.f90
index 51762c9262b2f0e7dec5b82bdcdfab64323438b9..59634c133c00cd136213de761e9c93a44d8d5a5f 100644
--- a/src/pbl_fo_turb.f90
+++ b/src/pbl_fo_turb.f90
@@ -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
diff --git a/src/pbl_solver.f90 b/src/pbl_solver.f90
index b5f00fde80a36d693fe95618de5471989fe06d1b..3313c4032739c053cdf30d81b0c0cf76978e55ea 100644
--- a/src/pbl_solver.f90
+++ b/src/pbl_solver.f90
@@ -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
diff --git a/src/pbl_turb_common.f90 b/src/pbl_turb_common.f90
index 437d4135e4f49020f88673bd2590017943960684..0d82b178c5db77d1620043dc2fafc833c7a9965f 100644
--- a/src/pbl_turb_common.f90
+++ b/src/pbl_turb_common.f90
@@ -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)
diff --git a/src/tests/gabls1.f90 b/src/tests/gabls1.f90
index f5c2bd1b9b8e8d72d000a3c3cdf3fb7967650b1e..76ffcc50bb2f89d107bf2e0ac652cd4af73ee32b 100644
--- a/src/tests/gabls1.f90
+++ b/src/tests/gabls1.f90
@@ -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