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

major update before merge

parent cdeb7427
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 418 additions and 83 deletions
config-examples/config-cbl-ex.txt
+ 5
− 4
View file @ 8fb231ed
......@@ -4,13 +4,14 @@ Vg = 0.0;
hbl_0 = 100.0;
Tgrad = 0.03;
fcoriolis = 1.39 * 0.0001;
tw_flux_bot = -0.35;
tw_flux_bot = 0.35;
time{
begin = 0.0;
end = 1.0 * 3600.0;
dt = 5.0;
out_dt = 360.0;
dt = 180.0;
out_dt = 1.0;
}
fluid {
......@@ -37,7 +38,7 @@ surface_model {
}
output{
dt = 0.2 *3600.0;
dt = 30.0; #0.2 *3600.0;
}
closure {
......
config-examples/config-gabls1-ex.txt
+ 1
− 1
View file @ 8fb231ed
......@@ -24,7 +24,7 @@ fluid {
grid {
type = "uniform";
nz = 64;
nz = 8;
h_bot = 0.0;
h_top = 800.0;
}
......
config-examples/config-gabls2-ex.txt
+ 3
− 3
View file @ 8fb231ed
......@@ -6,7 +6,7 @@ fcoriolis = 1.39 * 0.0001;
time{
begin = 0.0;
end = 60.0 * 3600.0;
dt = 20.0;
dt = 60.0;
out_dt = 360.0;
}
......@@ -21,7 +21,7 @@ fluid {
grid {
type = "uniform";
nz = 32;
nz = 24;
h_bot = 0.0;
h_top = 2000.0;
}
......@@ -34,7 +34,7 @@ surface_model {
}
output{
dt = 0.2 *3600.0;
dt = 60.0;
}
closure {
......
src/diag_pbl.f90
+ 6
− 4
View file @ 8fb231ed
......@@ -46,7 +46,7 @@ contains
write(*,*) 'here_hpbl'
dz_low = z(kl) - z_surf
hdyn = max(dz_low , 0.1E0 * ustar/cor)
hdyn = min(dz_low , 0.5E0 * ustar/cor)
write(*,*) 'hdyn', hdyn
kpblc = kl
kpbld = kl
......@@ -134,7 +134,7 @@ contains
real,intent(out)::hpbl
integer,intent(out)::nkpbl
real,parameter:: g = 9.81
real,parameter:: Ric = 0.25
real,parameter:: Ric = 0.5
real,parameter:: upper_bound = 16000.0 !upper boundary of PBLH
real::Rib
real du, dv, dtvirt
......@@ -146,7 +146,7 @@ contains
do k=kl-1,2,-1
if(z(k).lt.upper_bound)then
du = (u(k)-u(kl)) * (u(k)-u(kl)) + (v(k)-v(kl)) * (v(k)-v(kl))
du = (u(k)) * (u(k)) + (v(k)) * (v(k))
dtvirt = theta(k) - theta(kl)
if(du <= du2min) du = 0.01
......@@ -156,16 +156,18 @@ contains
kpbl=k
if(Rib.ge.Ric) then
!write(*,*) 'Rib', k, rib, du, dtvirt
kpbl=k
hpbl = z(k)
exit
end if
endif
enddo
do k=kl-1,2,-1
dtvirt = theta(k) - theta(kl)
if (dtvirt > 0 .and. kpbld == kl ) kpbld = k
end do
kpbl = amin0(kpbl, kpbld, kl-2)
kpbl = min(kpbl, kpbld+1, kl-2)
hpbl = z(kpbl)
hpbl = hpbl - zs
nkpbl = kpbl
......
src/fo_stability_functions.f90
+ 43
− 3
View file @ 8fb231ed
......@@ -21,14 +21,14 @@ module fo_stability_functions
!>
subroutine louis_stab_functions(fnriuv,fnritq,rinum,vdlm, &
& vdlh,zkh, rolim)
& vdlh,zkh, oc_flag)
implicit none
real fnriuv,fnritq
real, intent(in):: rinum
real, intent(in):: vdlm,vdlh, zkh
real scf, ucfm, ucfh
!real yb, yc, yd, ye, ricr
real, intent(in):: rolim
integer, intent(in):: oc_flag
if(rinum.ge.0.e0) then
scf = sqrt(1.e0 + yd*rinum)
......@@ -38,7 +38,7 @@ module fo_stability_functions
ucfm = 1.e0/(1.e0 + ye * (vdlm/zkh)**2 * sqrt(abs(rinum)))
ucfh = 1.e0/(1.e0 + ye * (vdlh/zkh)**2 * sqrt(abs(rinum)))
fnriuv = 1.e0 - 2.e0*yb * rinum * ucfm
if(rolim.gt.0.99.and.rolim.lt.1.5) then
if(oc_flag > 0) then
fnritq = 1.e0 - 3.0*3.e0*yb * rinum * ucfh !enhanced mixing for unstable stratification over land
else
fnritq = 1.e0 - 3.e0*yb * rinum * ucfh
......@@ -87,5 +87,45 @@ module fo_stability_functions
end if
end subroutine esau_vit_stab_functions
!>esau byrkjedal 2007 + viterbo 1999 unstable stratification
subroutine esau_stab_functions( fnriuv,fnritq,rinum, &
& vdlm, vdlh, zkh, dz, oc_flag)
implicit none
real fnriuv,fnritq,rinum, vdlm, vdlh, zkh, dz
integer oc_flag
real scf, ucfm, ucfh
real yb, yc, yd, ye, ricr
real gh
real, parameter:: p_m = 21.0
real, parameter:: q_m = 0.005
real, parameter:: p_h = 10.0
real, parameter:: q_h = 0.0012
real, parameter:: bb = 5.0
real, parameter:: cc = 5.0
yb = 5.e0
yc = 5.e0
yd = 5.e0
ye = yb*yc*sqrt(1.e0/3.e0)
ricr = 2.e0/(3.e0*yd)
if(rinum.ge.0.e0) then
scf = sqrt(1.e0 + yd*rinum)
fnriuv = 1.0/((1.0 + p_m * rinum)*(1.0 + p_m * rinum)) &
& + q_m * sqrt(rinum)
fnritq = 1.0/((1.0+p_h*rinum)*(1.0+p_h*rinum)*(1.0+p_h*rinum)) &
& + q_h
else
ucfm = 1.e0/(1.e0 + ye * (vdlm/zkh)**2 * sqrt(abs(rinum)))
ucfh = 1.e0/(1.e0 + ye * (vdlh/zkh)**2 * sqrt(abs(rinum)))
fnriuv = 1.e0 - 2.e0*yb * rinum * ucfm
if(oc_flag >= 1) then
fnritq = 1.e0 - 3.0*3.e0*yb * rinum * ucfh !enhanced mixing for unstable stratification over land
else
fnritq = 1.e0 - 3.e0*yb * rinum * ucfh
end if
end if
end subroutine esau_stab_functions
end module fo_stability_functions
\ No newline at end of file
src/pbl_dry_contrgradient.f90
+ 279
− 45
View file @ 8fb231ed
......@@ -16,9 +16,14 @@ module pbl_dry_contrgradient
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/3.0; !< vertical velocity in updraft Bernulli equation constant
real, parameter :: c_entr = 0.38 !< dry entrainment parameterization
real, parameter :: b1 = 1.8; !< vertical velocity in updraft Bernulli equation constant
real, parameter :: b2 = 3.5; !< vertical velocity in updraft Bernulli equation constant
real, parameter :: wstrmin = 1.0e-6 !< minimimal deardorff velocity
real, parameter :: lmin_entr = 1.0e6 !< maximum entrainement scale above CBL
real, parameter :: tau_dry_scale = 500.0 !< Time scale for entrainment
real, parameter :: a0 = 0.08 !< Updraft fraction area
real, parameter :: alph0 = 0.25 !< amount of surface dispersion to transfer to updraft
public cbl_allocate, cbl_deallocate
......@@ -64,14 +69,19 @@ module pbl_dry_contrgradient
integer k
real dz
cbl%entr(:) = 0.0
do k = grid%kmax, 2, -1
dz =(grid%z_cell(k-1) - grid%z_cell(k))
if (k <= bl%kpbl) then
dz =(grid%z_cell(k-1) - grid%z_cell(k)) + 1.0
if (k < bl%kpbl+2) 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));
cbl%entr(k) = c_entr * ( 1.0/(grid%z_cell(k) + dz) &
+ 1.0/(bl%hpbla_diag - grid%z_cell(k) + dz) &
+ 1.0/ lmin_entr)
end if
write(*,*) 'entr, k = ' , cbl%entr(k), k, grid%z_cell(k)
end do
cbl%entr(1) = 0
end subroutine cbl_get_entrainment
......@@ -86,51 +96,68 @@ module pbl_dry_contrgradient
type(stateBLDataType), intent(inout):: bl
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
real phys_beta, wstr, dz
real phys_beta, wstr, dz, dzentrp, dzentrm, temp
integer k, kmax
cbl%theta_up(:) = 0.0
cbl%thetav_up(:) = 0
cbl%qv_up = 0.0
cbl%ua_up = 0.0
cbl%va_up = 0.0
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)) &
write(*,*) 'bl%hs, ', (bl%surf%hs), phys_beta
if (bl%surf%hs > 0.0) then
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
wstr = max(wstrmin, wstr)
cbl%thetav_up0 = (bl%surf%hs/(bl%rho(kmax)*fluid%cp)) / wstr
cbl%theta_up0 = (bl%surf%hs/(bl%rho(kmax)*fluid%cp)) / wstr
cbl%qv_up0 = (bl%surf%es/bl%rho(kmax)) / wstr
else
cbl%thetav_up0 = 0.0
cbl%theta_up0 = 0.0
cbl%qv_up0 = 0.0
end if
cbl%ua_up0 = 0.0
cbl%va_up0 = 0.0
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
cbl%thetav_up(kmax) = cbl%thetav_up0 + bl%theta_v(kmax)
cbl%theta_up(kmax) = cbl%theta_up0 + bl%theta(kmax)
cbl%qv_up(kmax) = cbl%qv_up0 + bl%qv(kmax)
cbl%ua_up(kmax) = cbl%ua_up0 + bl%u(kmax)
cbl%va_up(kmax) = cbl%va_up0 + bl%v(kmax)
do k = kmax-1 ,1,-1
if (k <= bl%kpbl ) then
cbl%theta_up(k) = 0.0
cbl%thetav_up(k) = 0.0
cbl%qv_up(k) = 0
cbl%theta_up(k) = bl%theta(k)
cbl%thetav_up(k) = bl%theta_v(k)
cbl%qv_up(k) = bl%qv(k)
cbl%ua_up(k) = bl%u(k)
cbl%ua_up(k) = bl%u(k)
cbl%va_up(k) = bl%v(k)
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))
cbl%thetav_up(k) = cbl%thetav_up(k+1) - &
dz * cbl%entr(k+1) &
* (cbl%thetav_up(k+1) - bl%theta_v(k+1))
cbl%theta_up(k) = cbl%theta_up(k+1) - &
dz * cbl%entr(k+1) &
* (cbl%theta_up(k+1) - bl%theta(k+1))
write(*,*) 'theta_up, k = ' , cbl%theta_up(k) - bl%theta(k), k
cbl%qv_up(k) = cbl%qv_up(k+1) - &
dz * cbl%entr(k+1) &
* (cbl%qv_up(k+1) - bl%qv(k+1))
cbl%ua_up(k) = cbl%ua_up(k+1) + &
dz * cbl%entr(k+1) &
* (cbl%ua_up(k+1) - bl%u(k+1))
cbl%va_up(k) = cbl%va_up(k+1) - &
dz * cbl%entr(k+1) &
* (cbl%va_up(k+1) - bl%v(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
......@@ -204,7 +231,7 @@ module pbl_dry_contrgradient
!apply upwind
rhs_up(:) = 0
do k = kmax-1, 1, -1
if (abs(cbl%wu(k)) > 1.0e-8) then
if (abs(cbl%wu(k)) > 1.0e-7) 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))
......@@ -217,7 +244,7 @@ module pbl_dry_contrgradient
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)) * &
rhs_up(k) = 0.1 * (cbl%wu(k)) * &
(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
......@@ -225,11 +252,11 @@ module pbl_dry_contrgradient
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)) * &
rhs_up(k) = 0.1 * (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
......@@ -240,7 +267,7 @@ 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.00 * (cbl%wu(k)) * &
rhs_up(k) = 0.1 * (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
......@@ -251,4 +278,211 @@ module pbl_dry_contrgradient
end subroutine cbl_apply_cntrg
subroutine new_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
! local variables
real, save, allocatable :: rhs_up(:), rhs(:)
real, save, allocatable :: bouf(:)
real wu, phys_beta, umod, rhs_ws, a, b
real dz, dzentrp, dzentrm, temp,dth
real sigmaw, sigmath, sigmaqc
real wstr, ustr, thstar, thvstar, qstar
real z1, hpb, hpbmf
real entrint, entexp, entexpu
real wp, entw, wn2
real a0w0
integer kmax, kpbl, k
kmax = grid%kmax
kmax = grid%kmax
if (.not.(allocated(rhs_up))) then
allocate(rhs_up(grid%kmax), source=0.0)
end if
if (.not.(allocated(rhs))) then
allocate(rhs(grid%kmax), source=0.0)
end if
if (.not.(allocated(bouf))) then
allocate(bouf(grid%kmax), source=0.0)
end if
kpbl = bl%kpbl
cbl%theta_up(:) = 0.0
cbl%thetav_up(:) = 0
cbl%qv_up(:) = 0.0
cbl%ua_up(:) = 0.0
cbl%va_up(:) = 0.0
cbl%entr(:) = 0.0
cbl%wu(:) = 0.0
rhs_up(:) = 0.0
rhs(:) = 0.0
! Check if the flux is convective
if (bl%surf%hs > 0.0) then
! get surface properties
z1 = grid%z_cell(kmax)
hpb = bl%hpbla_diag
phys_beta = 0.5 * fluid%g /(bl%theta_v(kmax) + bl%theta_v(kmax-1))
wstr = (phys_beta * (bl%surf%hs/bl%rho(kmax)/fluid%cp) &
* bl%hpbla_diag)**(0.33333)
wstr = max(wstrmin, wstr)
! Initial condition for vertical velocity
umod = sqrt(bl%u(kmax)*bl%u(kmax) + bl%v(kmax)*bl%v(kmax))
ustr = sqrt(bl%surf%cm2u /bl%rho(kmax) * umod)
!cbl%wu(kmax) = 1.0 * ((ustr)**3.0 & !2.46
! + 0.4 * wstr**3.0)**(1.0/3.0) &
! * sqrt(1.0 - grid%z_cell(kmax) / bl%hpbla_diag)
sigmaw = 1.3 * (ustr**3.0 + 0.6 * wstr**3.0 * z1/hpb) * &
sqrt(1.0 - z1 / hpb)
thstar = (bl%surf%hs/(bl%rho(kmax))) / wstr
thvstar = (bl%surf%hs/(bl%rho(kmax))) / wstr
qstar = (bl%surf%es/bl%rho(kmax)) / wstr
cbl%wu(kmax) = 0.5 * sigmaw
! initial perturbation
cbl%thetav_up0 = (bl%surf%hs/(bl%rho(kmax))) / sigmaw
cbl%theta_up0 = (bl%surf%hs/(bl%rho(kmax))) / sigmaw
cbl%qv_up0 = (bl%surf%es/bl%rho(kmax)) / sigmaw
cbl%ua_up0 = 0.5 * (bl%u(kmax) + bl%u(kmax-1))
cbl%va_up0 = 0.5 * (bl%v(kmax) + bl%v(kmax-1))
cbl%thetav_up(kmax) = alph0 *cbl%thetav_up0 + bl%theta_v(kmax)
cbl%theta_up(kmax) = alph0 * cbl%theta_up0 + bl%theta(kmax)
cbl%qv_up(kmax) = alph0 * cbl%qv_up0 + bl%qv(kmax)
cbl%ua_up(kmax) = cbl%ua_up0 + bl%u(kmax)
cbl%va_up(kmax) = cbl%va_up0 + bl%v(kmax)
!Initial entrainment
!cbl%entr(kmax) = 1.0 / (cbl%wu(kmax) * tau_dry_scale) &
!+ c_entr / grid%z_cell(kmax)
dz = grid%z_cell(kmax) - grid%z_edge(kmax)
cbl%entr(kmax) = c_entr * ( 1.0/(grid%z_cell(kmax) + dz) &
+ 1.0/(bl%hpbla_diag - grid%z_cell(kmax) + dz))
bl%kpbld_mf = bl%kpbl
! get non-entraiened parcel height
!do k = kmax - 1, 2, -1
! dth = cbl%thetav_up(kmax) - bl%theta_v(k)
! if (bl%kpbld_mf == 0 .and. dth < 0.0 ) then
! bl%kpbld_mf = k
! end if
!end do
hpbmf = 0.5 * (grid%z_cell(bl%kpbld_mf) + &
bl%hpbla_diag)
write(*,*) 'hbp, hmf', grid%z_cell(bl%kpbld_mf), bl%hpbla_diag, bl%surf%hs
! get proper entrainment
do k = grid%kmax-1,1,-1
dz = (grid%z_cell(k) - grid%z_cell(k+1))
if (grid%z_cell(k) <= hpbmf) then
cbl%entr(k) = c_entr * ( 1.0/(grid%z_cell(k) + dz) &
+ 1.0/(hpbmf - grid%z_cell(k) + dz))
else
cbl%entr(k) = 0.0
end if
end do
!get updrafts
do k = grid%kmax-1,bl%kpbld_mf,-1
dz = (grid%z_cell(k) - grid%z_cell(k+1))
entexp = exp( - cbl%entr(k) * dz)
entexpu = exp( - cbl%entr(k) * dz / 3.0 )
cbl%thetav_up(k) = bl%theta_v(k) * (1.0 - entexp) + &
cbl%thetav_up(k+1) * entexp
cbl%theta_up(k) = bl%theta(k) * (1.0 - entexp) + &
cbl%theta_up(k+1) * entexp
cbl%qv_up(k) = bl%qv(k) * (1.0 - entexp) + &
cbl%qv_up(k+1) * entexp
cbl%ua_up(k) = bl%u(k) * (1.0 - entexpu) + &
cbl%ua_up(k+1) * entexpu
cbl%va_up(k) = bl%u(k) * (1.0 - entexpu) + &
cbl%va_up(k+1) * entexpu
bouf(k) = fluid%g * &
(0.5 * (cbl%thetav_up(k) + cbl%thetav_up(k+1)) / &
bl%theta_v(k) - 1.0)
enddo
!get vertical speed
do k = grid%kmax-1,bl%kpbld_mf,-1
dz = (grid%z_cell(k) - grid%z_cell(k+1))
wp = b1 * cbl%entr(k)
if (wp > 0.0) then
entw = exp( -2.0 * wp * dz )
wn2 = cbl%wu(k+1) * cbl%wu(k+1) * entw + &
b2* bouf(k)*(1.0-entw)/entw
else
wn2 = cbl%wu(k+1) * cbl%wu(k+1) + &
2.0 * b2 * bouf(k) * dz
end if
if ( wn2 > 0.0) then
cbl%wu(k) = sqrt(wn2)
else
cbl%wu(k) = 0.0
end if
end do
! apply upwind
rhs_up(:) = 0.0
do k = kmax-1, 1, -1
if (cbl%wu(k) > 0.0) then
rhs_up(k) = a0 * (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.0
end if
end do
bl%theta = bl%theta + dt * rhs_up
rhs_up(:) = 0.0
do k = kmax-1, 1, -1
if (cbl%wu(k) > 0.0) then
rhs_up(k) = a0 * (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.0
end if
end do
bl%qv(:) = bl%qv(:) + dt * rhs_up(:)
rhs_up(:) = 0.0
do k = kmax-1, 1, -1
if (cbl%wu(k) > 0.0) then
rhs_up(k) = a0 * (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.0
end if
end do
bl%u(:) = bl%u(:) + dt * rhs_up(:)
rhs_up(:) = 0.0
do k = kmax-1, 1, -1
if (cbl%wu(k) > 0.0) then
rhs_up(k) = a0 * (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.0
end if
end do
bl%v(:) = bl%v(:) + dt * rhs_up(:)
end if
end subroutine new_cntrg
end module pbl_dry_contrgradient
\ No newline at end of file
src/pbl_fo_turb.f90
+ 25
− 8
View file @ 8fb231ed
......@@ -22,6 +22,12 @@ module pbl_fo_turb
public get_fo_coeffs
public get_fo_coeffs_louis
public get_turb_length
type hpbl_options
integer:: hpbl_old = 1
integer:: hpbl_fixed = 2
integer:: hpbl_rib = 3
integer:: hpbl_heffner = 4
end type hpbl_options
! public :: default_stab_functions
! public :: esau_stab_functions
! private:: get_ph_efb
......@@ -70,14 +76,14 @@ module pbl_fo_turb
type(turbBLDataType), intent(inout):: turb
real fnriuv, fnritq,rho
real ocean_flag
integer ocean_flag
integer k, kmax
ocean_flag = real(bl%land_type)
ocean_flag = bl%land_type
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)
& turb%l_turb_h(k), grid%z_cell(k),ocean_flag)
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))
......@@ -105,14 +111,15 @@ module pbl_fo_turb
type(turbBLDataType), intent(inout):: turb
real fnriuv, fnritq,rho
real ocean_flag
integer ocean_flag
integer k, kmax
ocean_flag = real(bl%land_type)
ocean_flag = 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))
call esau_stab_functions( fnriuv,fnritq,turb%rig(k), &
turb%l_turb_m(k),turb%l_turb_h(k), grid%z_cell(k), &
grid%z_cell(k) - grid%z_cell(k+1), ocean_flag)
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))
......@@ -147,7 +154,12 @@ module pbl_fo_turb
integer k, kmax
kmax = grid%kmax
z_surf = grid%z_edge(kmax)
if (hbl_option == 1) then
if (hbl_option == 0) then
call get_hpbl(hpbla_diag, kpbl, bl%theta_v, grid%z_cell, &
z_surf, grid%kmax , bl%cor, turb%ustr)
bl%hpbla_diag = hpbla_diag
bl%kpbl = kpbl
elseif (hbl_option == 1) then
call get_hpbl(hpbla_diag, kpbl, bl%theta_v, grid%z_cell, &
z_surf, grid%kmax , bl%cor, turb%ustr)
bl%hpbla_diag = hpbla_diag
......@@ -157,6 +169,11 @@ module pbl_fo_turb
bl%u, bl%v, grid%kmax , grid%z_edge(kmax), du2min)
bl%hpbla_diag = hpbla_diag
bl%kpbl = kpbl
elseif (hbl_option == 3) then
call diag_pblh_rib(hpbla_diag,kpbl, bl%theta_v,grid%z_cell, &
bl%u, bl%v, grid%kmax , grid%z_edge(kmax), du2min)
bl%hpbla_diag = hpbla_diag
bl%kpbl = kpbl
else
write(*,*) 'wrong hbl diagnostics option'
return
......
src/pbl_solver.f90
+ 1
− 2
View file @ 8fb231ed
......@@ -146,7 +146,7 @@ module pbl_solver
ff(kmax) = bl%theta(kmax) &
+ (dt/grid%dze(kmax)) * bl%surf%hs / (bl%rho(kmax))
write(*,*) 'ff'!, ff!, aa,bb,cc
!write(*,*) 'ff'!, ff!, aa,bb,cc
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%theta)
......@@ -237,7 +237,6 @@ module pbl_solver
end do
! bottom boundary
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))
......
src/pbl_turb_common.f90
+ 1
− 1
View file @ 8fb231ed
......@@ -96,7 +96,7 @@ contains
call get_rig(turb, bl, fluid, grid, du2min, ricr)
call get_turb_length(turb, bl, fluid, grid, hbl_option)
write(*,*)'here5'
!write(*,*)'here5'
if (turb%cl_type == 1) then ! FOM closures
call get_fo_coeffs(turb, bl, fluid, grid)
......
src/scm_state_data.f90
+ 1
− 0
View file @ 8fb231ed
......@@ -29,6 +29,7 @@ module scm_state_data
real :: hpbla_diag !< diagnostic boundary layer height [m]
integer :: ktvdm !< top layer index for diffusion
integer :: kpbl !< boundary layer height index
integer :: kpbld_mf !< boundary layer height index for dry mf
integer :: land_type !< ocean/land flag
type(surfBLDataType) :: surf
......
src/tests/cbl_exp.f90
+ 45
− 6
View file @ 8fb231ed
......@@ -59,7 +59,7 @@ program cbl_exp
!io variables
real, dimension (5):: series_val
type (io_struct) :: series_f, scan_f
type (io_struct) :: series_f, scan_f, scan_cg_f
#ifdef USE_NETCDF
type(variable_metadata) :: meta_z, meta_t
type(variable_metadata) :: meta_temperature, meta_uwind, meta_vwind
......@@ -88,6 +88,10 @@ program cbl_exp
scan_f%fname='time_scan_cbl.dat'
call set_file(scan_f, scan_f%fname)
scan_cg_f%fname='time_scan_cgp_cbl.dat'
call set_file(scan_cg_f, scan_cg_f%fname)
#ifdef USE_NETCDF
call set_meta_nc(meta_z, meta_t, meta_temperature, meta_uwind, meta_vwind)
#endif
......@@ -289,7 +293,7 @@ program cbl_exp
write(*,*)'nt=0', nt
call to_file_1d_2var('temperature1.dat', grid%z_cell, bl%theta, bl_old%kmax)
do while (time_current < time_end)
call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
!call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
tsurf = bl_old%theta(kmax)
......@@ -299,7 +303,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
!meteo_cell%z0_t = z0/10
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
......@@ -307,7 +311,7 @@ program cbl_exp
bl%surf%cm2u = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U
taux = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U * bl_old%u(kmax)
tauy = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U * bl_old%v(kmax)
hflx = - bl%rho(kmax) &
hflx = bl%rho(kmax) &
* tw_flux0
bl%surf%hs = hflx
bl%surf%es = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U *meteo_cell%dQ
......@@ -316,14 +320,17 @@ program cbl_exp
!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)
call get_coeffs_general(turb, bl, fluid_params, 2, 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
! write(*,*) 'kdiffs', k, bl%vdctq(k), bl%vdcuv(k)
!end do
call new_cntrg(cbl, bl, fluid_params, grid, dt)
call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
!call cbl_apply_cntrg(cbl, bl, fluid_params, grid, dt)
call get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
!write(*,*) 'dttheta,', bl%theta- bl_old%theta
call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
!call to_file_1d_2var('temperature.dat', grid%z_cell, bl_old%U - bl%U, bl_old%kmax)
......@@ -352,6 +359,7 @@ program cbl_exp
call write_series(series_val, 5, series_f)
call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
call put_tscan_cg( time_current/3600.0, grid%z_cell,cbl, bl_old%kmax, scan_cg_f)
#ifdef USE_NETCDF
do k = 1, kl
ttemp(k,nt) = thold(k)
......@@ -459,6 +467,36 @@ subroutine put_tscan( time, z, bl, nl, f)
end subroutine put_tscan
subroutine put_tscan_cg( time, z, cbl, nl, f)
use parkinds, only : rf=>kind_rf, im => kind_im
use scm_io_default
use pbl_dry_contrgradient, only: pblContrGradDataType
implicit none
type(pblContrGradDataType), intent(in):: cbl
real(kind=rf), intent(in), dimension(nl):: z
real(kind=rf),intent(in) :: time
type (io_struct),intent(in) :: f
integer(kind=im), intent(in) :: nl
!local
real(kind=rf), dimension(5,nl)::stamp
integer k
! copy to stamp
do k=1,nl
stamp(1,k) = time
stamp(2,k) = z(k)
stamp(3,k) = cbl%wu(k)
stamp(4,k) = cbl%theta_up(k)
stamp(5,k) = cbl%entr(k)
end do
! call to write timestamp
call write_timescan(stamp,nl, 5, f)
end subroutine put_tscan_cg
subroutine get_surface_from_config(model, surface, z0,zh )
use iso_c_binding, only: C_NULL_CHAR
use config_parser
......@@ -609,3 +647,4 @@ subroutine set_meta_nc(z_meta, t_meta, theta_meta, uwind_meta, vwind_meta)
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
end subroutine set_meta_nc
#endif
src/tests/gabls1.f90
+ 1
− 1
View file @ 8fb231ed
......@@ -312,7 +312,7 @@ program gabls1
!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)
call get_coeffs_general(turb, bl, fluid_params, 0, 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/gabls2.f90
+ 7
− 5
View file @ 8fb231ed
......@@ -18,7 +18,7 @@ program gabls2
use diag_pbl
use pbl_grid
use sfx_data, only: meteoDataType, sfxDataType
use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
use sfx_esm, only: get_surface_fluxes_most => get_surface_fluxes, &
numericsType_most => numericsType
use config_parser
use config_utils
......@@ -254,7 +254,7 @@ program gabls2
meteo_cell%dQ = 0.0
meteo_cell%h = grid%z_cell(kmax)
meteo_cell%z0_m = z0
meteo_cell%z0_t= z0 / 10.0
!meteo_cell%z0_t= z0 / 10.0
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
......@@ -264,22 +264,24 @@ program gabls2
tauy = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U * bl_old%v(kmax)
hflx = - bl%rho(kmax) &
* sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U * meteo_cell%dT
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
!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)
call get_coeffs_general(turb, bl, fluid_params, 2, 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
! 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)
call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
call new_cntrg(cbl, bl, fluid_params, grid, dt)
!write(*,*) 'dttheta,', bl%theta - bl_old%theta
call scm_data_copy(bl_old,bl)
call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
!write(*,*) 'aftercoriolis,'
call get_wsub(w_sub, time_current, grid)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment