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

bug in tridiag

parent ed08e8fc
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
config-examples/config-gabls1-ex.txt
+ 2
− 1
View file @ ead9cb39
......@@ -4,6 +4,7 @@ Vg = 0.0;
hbl_0 = 100.0;
Tgrad = 0.01;
cooling_rate = -0.25 * 3600.0;
fcoriolis = 1.39 * 0.0001;
time{
begin = 0.0;
......@@ -23,7 +24,7 @@ fluid {
grid {
type = "uniform";
nz = 8;
nz = 16;
h_bot = 0.0;
h_top = 400.0;
}
......
src/diag_pbl.f90
+ 4
− 3
View file @ ead9cb39
......@@ -35,7 +35,7 @@ contains
implicit none
real, intent(in):: cor, ustar, z_surf
integer, intent(in):: kl
real, intent(in), dimension(KL):: theta_v, z
real, intent(in), dimension(*):: theta_v, z
real, intent(out):: hpbla
integer, intent(out):: kpbl
......@@ -44,10 +44,10 @@ contains
integer:: k, kpbld, kpblc
write(*,*) 'here_hpbl'
dz_low = z(kl) - z_surf
hdyn = AMIN1(dz_low , 0.5E0 * ustar/cor)
write(*,*) 'hdyn', hdyn
kpblc = kl
kpbld = kl
......@@ -59,6 +59,7 @@ contains
enddo
kpbl = MIN0(kpblc, kpbld, KL-2)
hpbla = z(kpbl) - z_surf
write(*,*) 'hpbla1, kpbl' , hpbla, kpbl
end subroutine get_hpbl
subroutine diag_pblh_inmcm(z,thetav,lat,zs,ustar,kl,hpbl)
......
src/pbl_fo_turb.f90
+ 22
− 11
View file @ ead9cb39
......@@ -74,19 +74,26 @@ module pbl_fo_turb
integer k, kmax
ocean_flag = real(bl%land_type)
do k = kmax, 1, -1
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)
turb%km(k) = fnriuv * turb%l_turb_m(k) * turb%l_turb_m(k) * turb%s2(k)
turb%kh(k) = fnritq * turb%l_turb_h(k) * turb%l_turb_h(k) * turb%s2(k)
!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
do k = kmax-1,1, -1
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
......@@ -102,13 +109,17 @@ module pbl_fo_turb
type(pblgridDataType), intent(in)::grid
type(turbBLDataType), intent(inout):: turb
integer hbl_option
real l_neutral, lm, lh
real l_neutral, lm, lh, hpbla_diag, z_surf
integer kpbl
integer k, kmax
kmax = grid%kmax
z_surf = grid%z_edge(kmax)
if (hbl_option ==1) then
call get_hpbl(bl%hpbla_diag, bl%kpbl, bl%theta_v, grid%z_cell, &
grid%z_edge(kmax), grid%kmax , bl%cor, turb%ustr)
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
else
write(*,*) 'wrong hbl diagnostics option'
return
......@@ -116,11 +127,11 @@ module pbl_fo_turb
do k = kmax, 1, -1
l_neutral = fluid%kappa * grid%z_cell(k)
if ( k <= bl%kpbl) then
if ( k < bl%kpbl+1) then
lm = 0.0
lh = 0.0
else
if (turb%rig(k) >= 0.0) then
if (turb%rig(k) > 0.0) then
lm = (1.0 - grid%z_cell(k) / bl%hpbla_diag)
lh = (1.0 - grid%z_cell(k) / bl%hpbla_diag)
else
......@@ -128,8 +139,8 @@ module pbl_fo_turb
lh = 1.0
endif
endif
turb%l_turb_m = l_neutral * lm
turb%l_turb_h = l_neutral * lh
turb%l_turb_m(k) = l_neutral * lm
turb%l_turb_h (k) = l_neutral * lh
end do
end subroutine get_turb_length
......
src/pbl_grid.f90
+ 2
− 1
View file @ ead9cb39
......@@ -98,9 +98,10 @@ module pbl_grid
grid%z_cell(nk) = grid%z_edge(nk) + 0.5 * dz
do k = nk-1, 1, -1
grid%z_edge(k) = grid%z_edge(k+1) + dz
grid%z_cell(k) = 0.5 * (grid%z_edge(k+1) + grid%z_edge(k))
grid%z_cell(k+1) = 0.5 * (grid%z_edge(k+1) + grid%z_edge(k))
end do
grid%z_edge(0) = h_top
grid%z_cell(1) = 0.5 * (grid%z_edge(1) + grid%z_edge(0))
grid%dzc(1:nk) = dz
grid%dze(1:nk) = dz
end subroutine set_pbl_grid_uniform
......
src/pbl_solver.f90
+ 81
− 13
View file @ ead9cb39
......@@ -40,70 +40,138 @@ module pbl_solver
integer, intent(in) :: ktvd, kl
real, intent(in), dimension(kl) :: aa, bb, cc, ff
real, intent(in), dimension(kl) :: prgna, prgnz
real, intent(out), dimension(kl) :: y
real, intent(inout), dimension(kl) :: y
integer :: k
real :: prgnb(kl)
write(*,*) 'here_diff3.1', bb(ktvd)
prgnb(ktvd) = ff(ktvd) / bb(ktvd)
write(*,*) 'here_diff3.1'
do k = ktvd+1, kl
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)
end do
end subroutine solve_tridiag
subroutine fill_tridiag(aa, bb, cc, rho, kdiff, kbltop, grid, dt)
subroutine fill_tridiag(aa, bb, cc, rho, kdiff, kbltop, cm2u, grid, dt)
use pbl_grid, only: pblgridDataType
implicit none
real, dimension(*), intent(in):: rho, kdiff
real, intent(in):: dt
real, intent(in):: dt, cm2u
integer, intent(in):: kbltop
type(pblgridDataType),intent(in):: grid
real, dimension(*), intent(out):: aa, bb, cc
real, dimension(*), intent(inout):: aa, bb, cc
real:: dtdz
integer:: k
!nulify before top boundary
aa(1:kbltop) = 0.0
bb(1:kbltop) = 0.0
cc(1:kbltop) = 0.0
aa(1:kbltop-1) = 0.0
bb(1:kbltop-1) = 0.0
cc(1:kbltop-1) = 0.0
!top boundary condition: flux = 0
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)
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)
end do
write(*,*) 'here_diff1.75'
!bottom boundary
k = grid%kmax
dtdz = dt / (grid%dzc(k))
aa(k) = (kdiff(k-1)/rho(k)) * dtdz / grid%dze(k-1)
bb(k) = 1.0 + aa(k)
bb(k) = 1.0 + aa(k) + dtdz * cm2u / rho(k)
cc(k) = 0.0
end subroutine fill_tridiag
subroutine solve_diffusion(bl, bl_old, turb, fluid, grid)
subroutine solve_diffusion(bl, bl_old, turb, 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(stateBLDataType), intent(out):: bl
type(stateBLDataType), intent(inout):: bl
type(stateBLDataType), intent(in):: bl_old
type(turbBLDataType), intent(in):: turb
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
real, intent(in):: dt
real, allocatable:: aa(:), bb(:), cc(:), ff(:)
real, allocatable:: prgna(:), prgnz(:)
integer k, integer, ktop, kmax
kmax = grid%kmax
write(*,*) 'here_diff'
if (.not.(allocated(aa))) then
allocate(aa(grid%kmax), source=0.0)
end if
if (.not.(allocated(bb))) then
allocate(bb(grid%kmax), source=0.0)
end if
if (.not.(allocated(cc))) then
allocate(cc(grid%kmax), source=0.0)
end if
if (.not.(allocated(ff))) then
allocate(ff(grid%kmax), source=0.0)
end if
if (.not.(allocated(prgna))) then
allocate(prgna(grid%kmax), source=0.0)
end if
if (.not.(allocated(prgnz))) then
allocate(prgnz(grid%kmax), source=0.0)
end if
write(*,*) 'here_diff1', bl%kpbl
ktop = bl%kpbl
do k = bl%kpbl-1,1,-1
if(bl%vdcuv(k) > 0.e0) then
ktop = k
end if
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'
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)
end do
ff(kmax) = bl%theta(kmax) + (dt/kmax) * bl%surf%hs /bl%rho(kmax)
write(*,*) 'here_diff3'
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%theta)
write(*,*) 'here_diff4'
do k = ktop,kmax-1
ff(k) = bl%qv(k)
end do
ff(kmax) = bl%qv(kmax) + (dt/kmax) * bl%surf%es /bl%rho(kmax)
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%qv)
write(*,*) 'here_diff4'
!velocity
call fill_tridiag(aa, bb, cc, bl%rho, bl%vdcuv, ktop, bl%surf%cm2u, grid, dt)
call factorize_tridiag(ktop, kmax, aa, bb, cc, prgna, prgnz)
do k = ktop,kmax-1
ff(k) = bl%u(k)
end do
ff(kmax) = bl%u(kmax)
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%u)
do k = ktop,kmax-1
ff(k) = bl%v(k)
end do
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
+ 3
− 5
View file @ ead9cb39
......@@ -43,7 +43,6 @@ contains
real du2
integer k
do k = 1,grid%kmax-1
write(*,*) 'k=', k
du2 = amax1(du2min, &
& (bl%u(k)-bl%u(k+1))**2 + (bl%v(k)-bl%v(k+1))**2)
turb%s2 = du2 / ((grid%z_cell(k) - grid%z_cell(k+1)) &
......@@ -66,7 +65,6 @@ contains
real dtvir, buoyp, du2, rinum
integer k
do k = 1,grid%kmax-1
write(*,*)'here rig'
du2 = amax1(du2min, &
& (bl%u(k)-bl%u(k+1))**2 + (bl%v(k)-bl%v(k+1))**2)
dtvir = bl%theta_v(k) - bl%theta_v(k+1)
......@@ -95,14 +93,14 @@ contains
call get_n2(turb, bl, fluid, grid)
write(*,*) 'here2'
call get_s2(turb, bl, fluid, grid, du2min)
write(*,*) 'here3'
call get_rig(turb, bl, fluid, grid, du2min, ricr)
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(*,*) 'here4'
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/scm_state_data.f90
+ 9
− 0
View file @ ead9cb39
......@@ -3,6 +3,11 @@
module scm_state_data
!< @brief abl state def.
implicit none
type surfBLDataType
real :: hs !< surface sensible heat flux W/m2
real :: es !< surface specific humidity flux kg/m2 s
real :: cm2u !< Cd * |U| for implicit momentum flux in tridiag solver
end type surfBLDataType
type stateBLDataType
real, allocatable:: u(:) !< u-componet velocity, [m/s]
real, allocatable:: v(:) !< v-componet velocity, [m/s]
......@@ -25,8 +30,10 @@ module scm_state_data
integer :: ktvdm !< top layer index for diffusion
integer :: kpbl !< boundary layer height index
integer :: land_type !< ocean/land flag
type(surfBLDataType) :: surf
end type stateBLDataType
public scm_data_allocate
public scm_data_deallocate
public scm_data_copy
......@@ -44,6 +51,7 @@ module scm_state_data
allocate(bl_data%qv(kmax))
allocate(bl_data%lwp(kmax))
allocate(bl_data%cloud_frac(kmax))
allocate(bl_data%rho(kmax))
allocate(bl_data%s_e(kmax))
allocate(bl_data%km(kmax))
allocate(bl_data%kh(kmax))
......@@ -62,6 +70,7 @@ module scm_state_data
deallocate(bl_data%qv)
deallocate(bl_data%lwp)
deallocate(bl_data%cloud_frac)
deallocate(bl_data%rho)
deallocate(bl_data%s_e)
deallocate(bl_data%km)
deallocate(bl_data%kh)
......
src/state_utils.f90
+ 23
− 2
View file @ ead9cb39
......@@ -184,10 +184,31 @@ module state_utils
real, intent(in):: p, th, qv
real t
t = (p / fluid%p0)**(fluid%R_d / fluid%cp)
get_density_theta = (p / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
t = th * (p / fluid%p0)**(fluid%R_d / fluid%cp)
get_density_theta = (p * 100.0 / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
end function get_density_theta
subroutine get_density_theta_vector( bl, fluid, grid)
use phys_fluid, only: fluidParamsDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
implicit none
type(fluidParamsDataType),intent(in):: fluid
type(pblgridDataType),intent(in):: grid
type(stateBLDataType),intent(inout)::bl
real t, p, qv
integer k
do k=grid%kmax, 1, -1
p = bl%p0 * grid%sigma(k)
t = bl%theta(k) * (p / fluid%p0)**(fluid%R_d / fluid%cp)
qv = bl%qv(k)
bl%rho(k) = (p * 100.0 / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%R_v/fluid%R_d)
end do
end subroutine get_density_theta_vector
subroutine get_theta_v( theta_v, theta, qv, fluid, kmax)
use phys_fluid, only: fluidParamsDataType
implicit none
......
src/tests/gabls1.f90
+ 23
− 8
View file @ ead9cb39
......@@ -13,7 +13,7 @@ program gabls1
use pbl_turb_common
use pbl_solver
use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z, &
get_sigma_from_z_theta, get_theta_v
get_sigma_from_z_theta, get_theta_v, get_density_theta_vector
use diag_pbl
use pbl_grid
use sfx_data, only: meteoDataType, sfxDataType
......@@ -198,6 +198,16 @@ program gabls1
stop
end if
end if
call c_config_is_varname("fcoriolis"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("fcoriolis"//C_NULL_CHAR, f_cor, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: fcoriolis '
ierr = 1 ! signal ERROR
stop
end if
end if
call set_fluid_default(fluid_params)
call get_fluid_params(fluid_params, status, ierr)
if (status == 0) then
......@@ -216,13 +226,14 @@ program gabls1
kmax = grid%kmax
time_current = time_begin
! hack to insure shir array is good
shir = 3.141592653589793238 * 72.3798 / 180.0
f_cor = 2.0 * fluid_params%omega * sin(shir)
!shir = 3.141592653589793238 * 72.3798 / 180.0
!f_cor = 2.0 * fluid_params%omega * sin(shir)
bl_old%cor= f_cor
bl%cor = f_cor
bl%land_type = 0.0
bl%p0=fluid_params%p0
write(*,*) 'F coriolis: = ', f_cor
!fs(1, 1) = 0.0
!setup surface
z0 = 0.01
zh = z0 / 10.0
......@@ -247,7 +258,7 @@ 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)
bl_old%theta_v(1:kmax) = bl_old%theta(1:kmax) * (1.0+ fluid_params%eps_v * bl_old%qv(1:kmax))
call get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
call scm_data_copy(bl,bl_old)
!Initialize turbulent closure
......@@ -280,15 +291,19 @@ program gabls1
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
cu = 1000.0 * sfx_cell%Cm**2
bl%surf%cm2u = sfx_cell%Cm**2 * meteo_cell%U
taux = 1.1 * sfx_cell%Cm**2 * meteo_cell%U * bl_old%u(kmax)
tauy = 1.1 * sfx_cell%Cm**2 * meteo_cell%U * bl_old%v(kmax)
hflx = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U * meteo_cell%dT
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
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 solv_diffusion(bl, bl_old, turb, fluid_params, 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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment