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

added gabls2 and cbl experiments

parent 12005e30
Branches main
Tags v.0.9
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1656 additions and 128 deletions
CMakeLists.txt
+ 42
− 0
View file @ e8953e1f
......@@ -100,18 +100,52 @@ if (WITH_TESTS)
target_link_libraries( gabls1 PRIVATE abl_lib )
target_link_libraries(gabls1 PRIVATE sfx_lib)
target_link_libraries(gabls1 PRIVATE config_parser_F config_parser_CXX)
#gabls2 experiment
add_executable(gabls2 src/tests/gabls2.f90
src/config-utils.f90
src/scm_io_default.f90
src/scm_sfx_data.f90)
target_include_directories(gabls2 PUBLIC ${CMAKE_BINARY_DIR}/modules/)
target_link_libraries( gabls2 PRIVATE abl_lib )
target_link_libraries(gabls2 PRIVATE sfx_lib)
target_link_libraries(gabls2 PRIVATE config_parser_F config_parser_CXX)
#cbl experiment
add_executable(cbl_exp src/tests/cbl_exp.f90
src/config-utils.f90
src/scm_io_default.f90
src/scm_sfx_data.f90)
target_include_directories(cbl_exp PUBLIC ${CMAKE_BINARY_DIR}/modules/)
target_link_libraries( cbl_exp PRIVATE abl_lib )
target_link_libraries(cbl_exp PRIVATE sfx_lib)
target_link_libraries(cbl_exp PRIVATE config_parser_F config_parser_CXX)
if(WITH_RRTMG)
target_link_libraries(gabls1 PRIVATE rrtm)
target_link_libraries(gabls2 PRIVATE rrtm)
target_link_libraries(cbl_exp PRIVATE rrtm)
endif()
target_compile_options(gabls1
PRIVATE $<$<COMPILE_LANGUAGE:Fortran>: -cpp>)
target_compile_options(gabls2
PRIVATE $<$<COMPILE_LANGUAGE:Fortran>: -cpp>)
target_compile_options(cbl_exp
PRIVATE $<$<COMPILE_LANGUAGE:Fortran>: -cpp>)
if ("${CMAKE_BUILD_TYPE}" STREQUAL "Debug")
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
target_compile_options(gabls1
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
target_compile_options(gabls2
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
target_compile_options(cbl_exp
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
target_compile_options(gabls1
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -O0 -fbacktrace -ffpe-trap=zero,overflow,underflow >)
target_compile_options(gabls2
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -O0 -fbacktrace -ffpe-trap=zero,overflow,underflow >)
target_compile_options(cbl_exp
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -O0 -fbacktrace -ffpe-trap=zero,overflow,underflow >)
endif()
endif()
......@@ -119,9 +153,17 @@ if (WITH_TESTS)
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
target_compile_options(gabls1
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
target_compile_options(gabls2
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
target_compile_options(cbl_exp
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
target_compile_options(gabls1
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -fbacktrace -ffpe-trap=zero,overflow,underflow >)
target_compile_options(gabls2
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -fbacktrace -ffpe-trap=zero,overflow,underflow >)
target_compile_options(cbl_exp
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -fbacktrace -ffpe-trap=zero,overflow,underflow >)
endif()
endif()
......
config-examples/config-cbl-ex.txt 0 → 100644
+ 46
− 0
View file @ e8953e1f
Ug = 7.5;
Vg = 0.0;
hbl_0 = 100.0;
Tgrad = 0.03;
fcoriolis = 1.39 * 0.0001;
tw_flux_bot = -0.35;
time{
begin = 0.0;
end = 1.0 * 3600.0;
dt = 5.0;
out_dt = 360.0;
}
fluid {
pref= 1013.4; #hPa
tref = 241.0; #K
g = 9.81; # m/s2
beta = g /tref;
kappa = 0.4; # von Karman constant
}
grid {
type = "uniform";
nz = 32;
h_bot = 0.0;
h_top = 1500.0;
}
surface_model {
id = "esm"; # optional: "esm" (default), "sheba", "most", "log"
type = "land";
z0_m = 0.1; # aerodynamic roughness [m]
z0_h = -1; # no prescribed value # -> using scheme assigned by surface type
}
output{
dt = 0.2 *3600.0;
}
closure {
type = "FOM";
name = "Louis"; # "Esau"; # "EFB";
}
\ No newline at end of file
config-examples/config-gabls1-ex.txt
+ 2
− 2
View file @ e8953e1f
......@@ -3,7 +3,7 @@ Ug = 8.0;
Vg = 0.0;
hbl_0 = 100.0;
Tgrad = 0.01;
cooling_rate = -0.25 * 3600.0;
cooling_rate = 0.25; # K/h
fcoriolis = 1.39 * 0.0001;
time{
......@@ -26,7 +26,7 @@ grid {
nz = 16;
h_bot = 0.0;
h_top = 400.0;
h_top = 800.0;
}
surface_model {
......
config-examples/config-gabls2-ex.txt 0 → 100644
+ 43
− 0
View file @ e8953e1f
Ug = 3.0;
Vg = -9.0;
fcoriolis = 1.39 * 0.0001;
time{
begin = 0.0;
end = 60.0 * 3600.0;
dt = 20.0;
out_dt = 360.0;
}
fluid {
pref= 1013.4; #hPa
tref = 283.15; #K
g = 9.81; # m/s2
beta = g /tref;
kappa = 0.4; # von Karman constant
}
grid {
type = "uniform";
nz = 32;
h_bot = 0.0;
h_top = 2000.0;
}
surface_model {
id = "esm"; # optional: "esm" (default), "sheba", "most", "log"
type = "land";
z0_m = 0.03; # aerodynamic roughness [m]
z0_h = -1; # no prescribed value # -> using scheme assigned by surface type
}
output{
dt = 0.2 *3600.0;
}
closure {
type = "FOM";
name = "Louis"; # "Esau"; # "EFB";
}
\ No newline at end of file
src/diag_pbl.f90
+ 54
− 52
View file @ e8953e1f
......@@ -22,7 +22,7 @@ contains
kpbld = kl
do k = kl-1,1,-1
dz_hdyn = dz_low - HDYN
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
......@@ -46,18 +46,18 @@ contains
write(*,*) 'here_hpbl'
dz_low = z(kl) - z_surf
hdyn = AMIN1(dz_low , 0.5E0 * ustar/cor)
hdyn = max(dz_low , 0.5E0 * ustar/cor)
write(*,*) 'hdyn', hdyn
kpblc = kl
kpbld = kl
do k = kl-1,1,-1
dz_hdyn = dz_low - HDYN
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
enddo
kpbl = MIN0(kpblc, kpbld, KL-2)
kpbl = MIN0(kpblc, kpbld, kl-2)
hpbla = z(kpbl) - z_surf
write(*,*) 'hpbla1, kpbl' , hpbla, kpbl
end subroutine get_hpbl
......@@ -114,59 +114,61 @@ contains
subroutine diag_pblh_rib(theta,z,u,kl,zs,hpbl,nkpbl)
! This subroutine calculates PBL depth according to (Troen and Mahrt 1986) as the lowest level where Rib>Ric
! Ric varies between 0.15 and 0.5
! Rib = (g/theta0)*(theta(z)-theta(s))*z/u(z)**2
!input:
! theta - array with (virtual) potential temperature
! z - array with heights above sea level
! u - array with wind speed
! zs - height of surface above sea level
! kl - number of vertical levels
! output:
! hpbl - PBL height above ground level
!!! IMPORTANT: In INMCM model vertical levels start form the top of the atmosphere, so arrays should be organized accordingly
implicit none
integer,intent(in)::kl
real,intent(in),dimension(kl)::z,theta,u
real,intent(in)::zs
real,intent(out)::hpbl
integer,intent(out)::nkpbl
real,parameter:: g = 9.81
real,parameter:: Ric = 0.25
real,parameter:: theta0 = 300.0
real,parameter:: upper_bound = 6000 !upper boundary of PBLH
real,dimension(kl)::Rib
real du
integer k,KPBL
KPBL=kl
do k=kl-1,1,-1
if(z(k).lt.upper_bound)then
du = u(k)-u(kl)
if(du.eq.0) du = 0.1
subroutine diag_pblh_rib(hpbl,nkpbl,theta,z,u,v,kl,zs, du2min)
! This subroutine calculates PBL depth according to (Troen and Mahrt 1986) as the lowest level where Rib>Ric
! Ric varies between 0.15 and 0.5
! Rib = (g/theta0)*(theta(z)-theta(s))*z/u(z)**2
!input:
! theta - array with (virtual) potential temperature
! z - array with heights above sea level
! u - array with wind speed
! zs - height of surface above sea level
! kl - number of vertical levels
! output:
! hpbl - PBL height above ground level
!!! IMPORTANT: In INMCM model vertical levels start form the top of the atmosphere, so arrays should be organized accordingly
implicit none
integer,intent(in)::kl
real,intent(in),dimension(kl)::z,theta,u,v
real,intent(in)::zs, du2min
real,intent(out)::hpbl
integer,intent(out)::nkpbl
real,parameter:: g = 9.81
real,parameter:: Ric = 0.25
real,parameter:: upper_bound = 16000.0 !upper boundary of PBLH
real::Rib
real du, dv, dtvirt
integer k,KPBL, kpbld, kpblc
Rib(k) = (g / theta0) * (theta(k) - theta(kl)) * (z(k) - z(kl)) / (du**2)
KPBL=kl
kpbld = kl
if(Rib(k).ge.Ric.and.KPBL.eq.kl)then
KPBL = k
if(KPBL.le.kl-2)then
hpbl = z(k) - (z(k)-z(k+1))*(Rib(k)-Ric)/(Rib(k)-Rib(k+1))
else
hpbl = z(k)
endif
endif
do k=kl-1,2,-1
if(z(k).lt.upper_bound)then
endif
enddo
du = (u(k)-u(kl)) * (u(k)-u(kl)) + (v(k)-v(kl)) * (v(k)-v(kl))
dtvirt = theta(k) - theta(kl)
if(du <= du2min) du = 0.01
if(KPBL.eq.kl) hpbl=z(kl)
Rib = (g * 2.0 / (theta(k) + theta(k+1))) &
* dtvirt * (z(k) - z(kl)) / du
hpbl = hpbl - zs
nkpbl = kl - KPBL + 1
kpbl=k
if(Rib.ge.Ric) then
!write(*,*) 'Rib', k, rib, du, dtvirt
hpbl = z(k)
exit
end if
endif
enddo
do k=kl-1,2,-1
if (dtvirt > 0 .and. kpbld == kl ) kpbld = k
end do
kpbl = amin0(kpbl, kpbld, kl-2)
hpbl = z(kpbl)
hpbl = hpbl - zs
nkpbl = kpbl
end subroutine diag_pblh_rib
......
src/pbl_dry_contrgradient.f90
+ 20
− 19
View file @ e8953e1f
......@@ -20,12 +20,13 @@ module pbl_dry_contrgradient
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 cbl_allocate, cbl_deallocate
public get_entrainment, get_up, get_wu, apply_cntrg
public cbl_get_entrainment, cbl_get_up
public cbl_get_wu, cbl_apply_cntrg
contains
subroutine allocate_cbl(cbl, kmax)
subroutine cbl_allocate(cbl, kmax)
integer, intent(in):: kmax
type(pblContrGradDataType), intent(out):: cbl
......@@ -36,9 +37,9 @@ module pbl_dry_contrgradient
allocate(cbl%ua_up(kmax))
allocate(cbl%va_up(kmax))
allocate(cbl%wu(kmax))
end subroutine allocate_cbl
end subroutine cbl_allocate
subroutine deallocate_cbl(cbl)
subroutine cbl_deallocate(cbl)
type(pblContrGradDataType), intent(inout):: cbl
deallocate(cbl%entr)
......@@ -48,8 +49,8 @@ module pbl_dry_contrgradient
deallocate(cbl%ua_up)
deallocate(cbl%va_up)
deallocate(cbl%wu)
end subroutine deallocate_cbl
subroutine get_entrainment(cbl, bl, fluid, grid)
end subroutine cbl_deallocate
subroutine cbl_get_entrainment(cbl, bl, fluid, grid)
use scm_state_data, only : stateBLDataType
use pbl_turb_data, only : turbBLDataType
use phys_fluid, only: fluidParamsDataType
......@@ -73,9 +74,9 @@ module pbl_dry_contrgradient
end if
end do
cbl%entr(1) = 0
end subroutine get_entrainment
end subroutine cbl_get_entrainment
subroutine get_up(cbl, bl, fluid, grid)
subroutine cbl_get_up(cbl, bl, fluid, grid)
use scm_state_data, only : stateBLDataType
use pbl_turb_data, only : turbBLDataType
use phys_fluid, only: fluidParamsDataType
......@@ -131,9 +132,9 @@ module pbl_dry_contrgradient
cbl%va_up(k) = cbl%va_up(k) + bl%v(k)
end do
end subroutine get_up
end subroutine cbl_get_up
subroutine get_wu(cbl, bl, fluid, grid)
subroutine cbl_get_wu(cbl, bl, fluid, grid)
use scm_state_data, only : stateBLDataType
use pbl_turb_data, only : turbBLDataType
use phys_fluid, only: fluidParamsDataType
......@@ -174,9 +175,9 @@ module pbl_dry_contrgradient
cbl%wu(k) = 0.0;
endif
enddo
end subroutine get_wu
end subroutine cbl_get_wu
subroutine apply_cntrg(cbl, bl, fluid, grid, dt)
subroutine cbl_apply_cntrg(cbl, bl, fluid, grid, dt)
use scm_state_data, only : stateBLDataType
use pbl_turb_data, only : turbBLDataType
use phys_fluid, only: fluidParamsDataType
......@@ -197,14 +198,14 @@ module pbl_dry_contrgradient
a0w0 = 0.15 * cbl%wu(kpbl)
call get_entrainment(cbl, bl, fluid, grid)
call get_up(cbl, bl, fluid, grid)
call get_wu(cbl, bl, fluid, grid)
call cbl_get_entrainment(cbl, bl, fluid, grid)
call cbl_get_up(cbl, bl, fluid, grid)
call cbl_get_wu(cbl, bl, fluid, grid)
!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)) * &
rhs_up(k) = 0.05 * (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
......@@ -216,7 +217,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.1 * (cbl%wu(k)) * &
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))
else
......@@ -224,6 +225,6 @@ module pbl_dry_contrgradient
end if
end do
bl%qv(:) = bl%qv(:) + dt * rhs_up(:)
end subroutine apply_cntrg
end subroutine cbl_apply_cntrg
end module pbl_dry_contrgradient
\ No newline at end of file
src/pbl_fo_turb.f90
+ 10
− 5
View file @ e8953e1f
......@@ -18,7 +18,7 @@ module pbl_fo_turb
! cormin is a minimal value of coriolis parameter to calculate the
! pbl extention
real, parameter:: du2min = 1.e-02
real, parameter:: cormin = 5.e-05
real, parameter:: cormin = 5.e-05 !include in hpbl diag
public get_fo_coeffs
public get_fo_coeffs_louis
public get_turb_length
......@@ -78,12 +78,12 @@ module pbl_fo_turb
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)
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%vdcuv(grid%kmax) = 0.0
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)
......@@ -110,11 +110,16 @@ module pbl_fo_turb
integer k, kmax
kmax = grid%kmax
z_surf = grid%z_edge(kmax)
if (hbl_option ==1) then
if (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
bl%kpbl = kpbl
elseif (hbl_option == 2) 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
+ 125
− 34
View file @ e8953e1f
......@@ -10,6 +10,7 @@ module pbl_solver
public solve_tridiag
public fill_tridiag
public solve_diffusion
public apply_subsidence
contains
subroutine factorize_tridiag(ktvd, kl, aa, bb, cc, prgna, prgnz)
implicit none
......@@ -44,11 +45,9 @@ module pbl_solver
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
y(kl) = prgnb(kl)
......@@ -71,31 +70,27 @@ module pbl_solver
integer:: k
!nulify before top boundary
aa(1:kbltop-1) = 0.0
bb(1:kbltop-1) = 0.0
cc(1:kbltop-1) = 0.0
aa(1:grid%kmax) = 0.0
bb(1:grid%kmax) = 0.0
cc(1:grid%kmax) = 0.0
!top boundary condition: flux = 0
!write(*,*) 'here_diff1.25', kbltop
k = kbltop
dtdz = dt / (grid%dzc(k))
dtdz = dt / (grid%z_edge(k) - grid%z_edge(k-1) )
aa(k) = 0
cc(k) = (kdiff(k)/rho(k)) * dtdz / grid%dze(k)
cc(k) = (kdiff(k)/rho(k)) * dtdz / (grid%z_cell(k+1) - grid%z_cell(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)
dtdz = dt / (grid%z_edge(k) - grid%z_edge(k-1) )
aa(k) = (kdiff(k - 1)/rho(k)) * dtdz / (grid%z_cell(k) - grid%z_cell(k-1))
cc(k) = (kdiff(k)/rho(k)) * dtdz /(grid%z_cell(k+1) - grid%z_cell(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) + dtdz * cm2u / rho(k)
dtdz = dt / (grid%z_edge(k) - grid%z_edge(k-1) )
aa(k) = (kdiff(k-1)/rho(k)) * dtdz / (grid%z_cell(k) - grid%z_cell(k-1))
bb(k) = 1.0 + aa(k) - dtdz * cm2u / rho(k)
cc(k) = 0.0
end subroutine fill_tridiag
......@@ -116,7 +111,6 @@ module pbl_solver
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
......@@ -135,45 +129,142 @@ module pbl_solver
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
!this loop is generally not needed
!do k = bl%kpbl-1,1,-1
! if(bl%vdcuv(k) > 0.e0) then
! ktop = k
! end if
!enddo
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'
ff(kmax) = bl%theta(kmax) &
+ (dt/grid%dze(kmax)) * bl%surf%hs / (bl%rho(kmax))
write(*,*) 'ff'!, ff!, aa,bb,cc
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)
ff(kmax) = bl%qv(kmax) + (dt/grid%dze(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
subroutine apply_subsidence(bl, w_s, 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(inout):: bl
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
real, intent(in):: dt
real, dimension(*), intent(in):: w_s
real, save, allocatable :: subs_u(:), subs_v(:), &
subs_theta(:), subs_qv(:)
real:: dtdz, rhom, rhop
integer k, kmax
kmax = grid%kmax
if (.not.(allocated(subs_u))) then
allocate(subs_u(grid%kmax), source=0.0)
end if
if (.not.(allocated(subs_v))) then
allocate(subs_v(grid%kmax), source=0.0)
end if
if (.not.(allocated(subs_theta))) then
allocate(subs_theta(grid%kmax), source=0.0)
end if
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))
rhom = 0.5*(bl%rho(k) + bl%rho(k-1))
subs_u(k) = (0.5/bl%rho(k)) &
* (rhop * w_s(k) * (bl%u(k)-bl%u(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1)) &
+ &
rhom * w_s(k-1) * (bl%u(k)-bl%u(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1)))
subs_v(k) = (0.5/bl%rho(k)) &
* (rhop * w_s(k) * (bl%v(k)-bl%v(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1)) &
+ &
rhom * w_s(k-1) * (bl%v(k)-bl%v(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1)))
subs_theta(k) = (0.5/bl%rho(k)) &
* (rhop * w_s(k) * (bl%theta(k)-bl%theta(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1)) &
+ &
rhom * w_s(k-1) * (bl%theta(k)-bl%theta(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1)))
subs_qv(k) = (0.5/bl%rho(k)) &
* (rhop * w_s(k) * (bl%qv(k)-bl%qv(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1)) &
+ &
rhom * w_s(k-1) * (bl%qv(k)-bl%qv(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1)))
end do
! bottom boundary
write(*,*) 'ws', grid%z_cell(kmax), grid%z_cell(kmax-1)
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))
subs_v(k) = w_s(k-1) * (bl%v(k)-bl%v(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1))
subs_theta(k) = w_s(k-1) * (bl%theta(k)-bl%theta(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1))
subs_qv(k) = w_s(k-1) * (bl%qv(k)-bl%qv(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1))
! top boundary
k = 1
subs_u(k) = w_s(k+1) * (bl%u(k)-bl%u(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1))
subs_v(k) = w_s(k+1) * (bl%v(k)-bl%v(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1))
subs_theta(k) = w_s(k+1) * (bl%theta(k)-bl%theta(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1))
subs_qv(k) = w_s(k+1) * (bl%qv(k)-bl%qv(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1))
!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
end do
end subroutine apply_subsidence
end module pbl_solver
\ No newline at end of file
src/pbl_turb_common.f90
+ 2
− 1
View file @ e8953e1f
......@@ -94,8 +94,9 @@ contains
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)
call get_turb_length(turb, bl, fluid, grid, hbl_option)
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 @ e8953e1f
......@@ -93,6 +93,7 @@ module scm_state_data
bl%qv = bl_old%qv
bl%lwp = bl_old%lwp
bl%cloud_frac= bl_old%cloud_frac
bl%rho= bl_old%rho
bl%s_e = bl_old%s_e
bl%km = bl_old%km
bl%kh = bl_old%kh
......
src/tests/cbl_exp.f90 0 → 100644
+ 610
− 0
View file @ e8953e1f
! Created by Andrey Debolskiy on 12.10.2024.
program cbl_exp
use ISO_C_BINDING, only: C_NULL_CHAR
use parkinds, only: rf=>kind_rf, im=>kind_im
use phys_fluid
use scm_io_default
#ifdef USE_NETCDF
use io, only: variable_metadata, write_2d_real_nc, write_1d_real_nc
#endif
use scm_sfx_data, only: taux, tauy, umst, hflx, qflx, cu
use scm_state_data
use pbl_turb_data
use pbl_turb_common
use pbl_dry_contrgradient
use pbl_solver
use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z, &
get_sigma_from_z_theta, get_theta_v, get_density_theta_vector
use diag_pbl
use pbl_grid
use sfx_data, only: meteoDataType, sfxDataType
use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
numericsType_most => numericsType
use config_parser
use config_utils
implicit none
!local varaibles
real time_begin, time_end, time_current
real dt, out_dt, output_dt
integer status, ierr
logical bstatus
character(len = 160) ::fname
character(len = 128) :: fname_nc
character(len = 128) ::fname_config = '../config-examples/config-cbl-ex.txt'
real, allocatable :: ttemp(:,:), utemp(:,:),vtemp(:,:), ttime(:)
real(kind=rf):: ug, vg
real(kind=rf):: z0, zh, f_cor
real(kind=rf):: tsurf
real(kind=rf):: shir
real(kind=rf):: hbl_0
real(kind=rf):: tgrad
real(kind=rf):: tw_flux0
type(fluidParamsDataType) fluid_params
type(pblgridDataType) grid
type(stateBLDataType):: bl, bl_old;
type(turbBLDataType):: turb
type(pblContrGradDataType):: cbl;
type(meteoDataType) :: meteo_cell
type(sfxDataType) :: sfx_cell
type(numericsType_most) :: numerics_sfx
!diagnostic variables
real hpbl
integer nkpbl
!io variables
real, dimension (5):: series_val
type (io_struct) :: series_f, scan_f
#ifdef USE_NETCDF
type(variable_metadata) :: meta_z, meta_t
type(variable_metadata) :: meta_temperature, meta_uwind, meta_vwind
#endif
integer k, nt, kmax
! init experiment params
numerics_sfx%maxiters_charnock = 10
ug = 8.00000
vg = 0.0
dt = 1.0
! call set_fluid_default(fluid_params)
! fluid_params%tref = 265.0
! fluid_params%pref = 1013.2500
! fluid_params%beta = fluid_params%g / fluid_params%tref
! fluid_params%kappa = 0.4
! fluid_params%p0= fluid_params%pref
!set output filenames
fname = 'test3.dat'
fname_nc = 'cbl.nc'
series_f%fname = 'phys_cbl.dat'
call set_file(series_f, series_f%fname)
scan_f%fname='time_scan_cbl.dat'
call set_file(scan_f, scan_f%fname)
#ifdef USE_NETCDF
call set_meta_nc(meta_z, meta_t, meta_temperature, meta_uwind, meta_vwind)
#endif
call init_config(fname_config,status, ierr)
if (status == 0) then
write(*, *) ' FAILURE! > could not initialize config file: ', fname_config
ierr = 1 ! signal ERROR
stop
end if
call c_config_is_varname("time.begin"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("time.begin"//C_NULL_CHAR, time_begin, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: time.begin '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("time.end"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("time.end"//C_NULL_CHAR, time_end, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: time.end '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("time.dt"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("time.dt"//C_NULL_CHAR, dt, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: time.dt '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("time.out_dt"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("output.dt"//C_NULL_CHAR, output_dt, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: output.dt '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("Ug"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("Ug"//C_NULL_CHAR, ug, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: Ug '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("Vg"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("Vg"//C_NULL_CHAR, vg, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: Vg '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("hbl_0"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("hbl_0"//C_NULL_CHAR, hbl_0, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: hbl_0 '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("Tgrad"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("Tgrad"//C_NULL_CHAR, tgrad, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: Tgrad '
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
write(*, *) ' FAILURE! > could not set: Fluid params '
ierr = 1 ! signal ERROR
stop
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 c_config_is_varname("tw_flux_bot"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("tw_flux_bot"//C_NULL_CHAR, tw_flux0, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: tw_flux_bot '
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
write(*, *) ' FAILURE! > could not set: Fluid params '
ierr = 1 ! signal ERROR
stop
end if
!z coordinate
call get_grid_params(grid, status, ierr)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: grid '
ierr = 1 ! signal ERROR
stop
end if
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)
bl_old%cor= f_cor
bl%cor = f_cor
bl%land_type = 0.0
bl%p0=fluid_params%p0
!setup surface
z0 = 0.01
zh = z0 / 10.0
! init blData states
call scm_data_allocate(bl, kmax)
call scm_data_allocate(bl_old, kmax)
call init_state(bl, ug, vg, fluid_params%tref)
call scm_data_copy(bl_old,bl)
call cbl_allocate(cbl, kmax)
DO k = 1, kmax
if (grid%z_cell(k) > 100.0) then
bl%theta(k) = fluid_params%tref + tgrad * (grid%z_cell(k) - 100.0)
bl%u(k) = ug
bl%v(k) = vg
else
bl%theta(k) = fluid_params%tref
!bl%u(k) = ug * 0.5
!bl%v(k) = vg * 0.5
end if
!write(*,*) k, sig(k)* fluid_params%pref* 100.0, grid%z_cell(k)
END DO
!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_old,bl)
!Initialize turbulent closure
call turb_data_allocate(turb, kmax)
! getclosurefromconfig
call get_closure_params(turb, status, ierr)
!io setup
nt = floor((time_end - time_begin)/output_dt)
allocate(ttime(nt))
allocate(ttemp(kmax,nt),utemp(kmax,nt),vtemp(kmax,nt))
#ifdef USE_NETCDF
call write_1d_real_nc(grid%z_cell,fname_nc,meta_z)
#endif
out_dt = 0.0
nt=0
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)
tsurf = bl_old%theta(kmax)
meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%u(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
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
call get_density_theta_vector( bl, fluid_params, grid)
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) &
* tw_flux0
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, 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)
!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)
call scm_data_copy(bl_old,bl)
time_current = time_current + dt
if (time_current >=out_dt) then
nt = nt+1
ttime(nt) = time_current/3600.0
out_dt = out_dt + output_dt
write(*,*)'nt= ', nt
call theta2ta(bl_old%temp, bl_old%theta, &
fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call diag_pblh_inmcm(grid%z_cell,bl_old%theta,shir,0.0,umst,bl_old%kmax,hpbl)
!create output
series_val(1) = time_current/3600.0
series_val(2) = hflx
series_val(3) = turb%ustr
series_val(4) = hpbl
series_val(5) = 0.0
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)
#ifdef USE_NETCDF
do k = 1, kl
ttemp(k,nt) = thold(k)
utemp(k,nt) = uold(k)
vtemp(k,nt) = vold(k)
end do
#endif
!write(*,*) ccld * 10000.0
end if
end do
#ifdef USE_NETCDF
write(*,*)'nt= ', nt
call write_1d_real_nc(ttime,fname_nc,meta_t)
write(*,*)'here'
call write_2d_real_nc(ttemp,fname_nc,meta_temperature)
call write_2d_real_nc(utemp,fname_nc,meta_uwind)
call write_2d_real_nc(vtemp,fname_nc,meta_vwind)
#endif
call ta2theta(bl_old%theta, bl_old%temp, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call to_file_1d_3var(fname, grid%z_cell, bl_old%u, bl_old%v, bl_old%kmax)
!call to_file_1d_3var('coeffs.dat',grid%z_cell, bl_old%km, bl_old%kh, bl_old%kmax)
call to_file_1d_2var('temperature.dat', grid%z_cell, bl%theta, bl_old%kmax)
call close_file(series_f)
call close_file(scan_f)
call deallocate_pbl_grid(grid)
call scm_data_deallocate(bl)
call scm_data_deallocate(bl_old)
call pbl_data_deallocate(turb)
call cbl_deallocate(cbl)
end program cbl_exp
subroutine init_state(bl, ug_,vg_,tref)
use parkinds, only: rf=>kind_rf, im=>kind_im
use scm_state_data, only:stateBLDataType
implicit none
real(kind=rf):: ug_,vg_,tref
type(stateBLDataType), intent(inout) :: bl
integer kmax
kmax= bl%kmax
write(*,*) 'blkmax=', kmax
bl%u(1:bl%kmax) = ug_
bl%v(1:bl%kmax) = vg_
bl%temp(1:bl%kmax) = tref
bl%theta(1:bl%kmax) = tref
bl%qv(1:bl%kmax) = 0.0
bl%lwp(1:bl%kmax) = 0.0
bl%cloud_frac(1:bl%kmax) = 0.0
bl%s_e(1:bl%kmax) = 0.0
bl%km(1:bl%kmax) = 0.0
bl%kh(1:bl%kmax) = 0.0
bl%vdcuv(1:bl%kmax) = 0.0
bl%vdctq(1:bl%kmax) = 0.0
end subroutine init_state
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(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)
end do
end subroutine add_coriolis
subroutine put_tscan( time, z, bl, nl, f)
use parkinds, only : rf=>kind_rf, im => kind_im
use scm_io_default
use scm_state_data, only:stateBLDataType
implicit none
type(stateBLDataType), intent(in):: bl
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) = bl%theta(k)
stamp(4,k) = bl%u(k)
stamp(5,k) = bl%v(k)
end do
! call to write timestamp
call write_timescan(stamp,nl, 5, f)
end subroutine put_tscan
subroutine get_surface_from_config(model, surface, z0,zh )
use iso_c_binding, only: C_NULL_CHAR
use config_parser
use sfx_common, only: char_array2str
use sfx_config
use sfx_surface, only: get_surface_id
use parkinds, only: rf=>kind_rf
implicit none
character, allocatable :: config_field(:)
integer,intent(out) :: model, surface
real(kind=rf), intent(out) :: z0,zh
integer status, ierr
call c_config_is_varname("model.id"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_string("model.id"//C_NULL_CHAR, config_field, status)
if (status == 0) then
ierr = 1 ! signal ERROR
stop
end if
model = get_model_id(char_array2str(config_field))
if (model == -1) then
write(*, *) ' FAILURE! > unknown model [key]: ', trim(char_array2str(config_field))
ierr = 1 ! signal ERROR
return
end if
end if
call c_config_is_varname("surface.type"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_string("surface.type"//C_NULL_CHAR, config_field, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
surface = get_surface_id(char_array2str(config_field))
if (surface == -1) then
write(*, *) ' FAILURE! > unknown surface [key]: ', trim(char_array2str(config_field))
ierr = 1 ! signal ERROR
return
end if
endif
end subroutine get_surface_from_config
#ifdef USE_NETCDF
subroutine set_meta_nc(z_meta, t_meta, theta_meta, uwind_meta, vwind_meta)
use io, only: variable_metadata
type(variable_metadata), intent(inout):: z_meta, t_meta, theta_meta, uwind_meta, vwind_meta
z_meta = variable_metadata( &
name = 'Z', &
dim_names = [character(len=32) :: 'Z', '', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Height', &
'm', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
t_meta = variable_metadata( &
name = 'Time', &
dim_names = [character(len=32) :: 'Time', '', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Time', &
'hours since 00-00-00', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
theta_meta = variable_metadata( &
name = 'th', &
dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Potential Temperature', &
'K', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
uwind_meta = variable_metadata( &
name = 'ua', &
dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Longitude wind component', &
'm/s', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
vwind_meta = variable_metadata( &
name = 'va', &
dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Latitude wind component', &
'm/s', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
end subroutine set_meta_nc
#endif
\ No newline at end of file
src/tests/gabls1.f90
+ 30
− 15
View file @ e8953e1f
......@@ -68,7 +68,7 @@ program gabls1
ug = 8.00000
vg = 0.0
dt = 1.0;
dt = 1.0
! call set_fluid_default(fluid_params)
! fluid_params%tref = 265.0
! fluid_params%pref = 1013.2500
......@@ -247,11 +247,15 @@ program gabls1
DO K = 1, kmax
DO k = 1, kmax
if (grid%z_cell(k) > 100.0) then
bl_old%theta(k) = fluid_params%tref + tgrad * (grid%z_cell(k) - 100.0)
bl%theta(k) = fluid_params%tref + tgrad * (grid%z_cell(k) - 100.0)
bl%u(k) = ug
bl%v(k) = vg
else
bl_old%theta(k) = fluid_params%tref
bl%theta(k) = fluid_params%tref
!bl%u(k) = ug * 0.5
!bl%v(k) = vg * 0.5
end if
!write(*,*) k, sig(k)* fluid_params%pref* 100.0, grid%z_cell(k)
END DO
......@@ -280,7 +284,10 @@ program gabls1
out_dt = 0.0
nt=0
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)
tsurf = fluid_params%tref - cooling_rate * (time_current - time_begin) / 3600.0;
meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%u(kmax)**2)
......@@ -292,21 +299,29 @@ program gabls1
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
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
call get_density_theta_vector( bl, fluid_params, grid)
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) &
* 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
umst = turb%ustr
call get_density_theta_vector( bl, fluid_params, grid)
!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
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)
!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)
call scm_data_copy(bl_old,bl)
......@@ -318,7 +333,7 @@ program gabls1
write(*,*)'nt= ', nt
call ta2theta( bl_old%theta, bl_old%temp, &
call theta2ta(bl_old%temp, bl_old%theta, &
fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call diag_pblh_inmcm(grid%z_cell,bl_old%theta,shir,0.0,umst,bl_old%kmax,hpbl)
......@@ -353,8 +368,8 @@ program gabls1
call ta2theta(bl_old%theta, bl_old%temp, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call to_file_1d_3var(fname, grid%z_cell, bl_old%u, bl_old%v, bl_old%kmax)
call to_file_1d_3var('coeffs.dat',grid%z_cell, bl_old%km, bl_old%kh, bl_old%kmax)
call to_file_1d_2var('temperature.dat', grid%z_cell, bl_old%theta, bl_old%kmax)
!call to_file_1d_3var('coeffs.dat',grid%z_cell, bl_old%km, bl_old%kh, bl_old%kmax)
call to_file_1d_2var('temperature.dat', grid%z_cell, bl%theta, bl_old%kmax)
call close_file(series_f)
call close_file(scan_f)
......@@ -402,8 +417,8 @@ subroutine add_coriolis(bl, bl_old, ug, vg, f, dt, grid)
integer k
do k = 1, grid%kmax
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)
bl%v(k) = bl%v(k) - dt * f * (bl%u(k) - ug)
bl%u(k) = bl%u(k) + dt * f * (bl%v(k) - vg)
end do
end subroutine add_coriolis
......
src/tests/gabls2.f90 0 → 100644
+ 671
− 0
View file @ e8953e1f
! Created by Andrey Debolskiy on 12.10.2024.
program gabls2
use ISO_C_BINDING, only: C_NULL_CHAR
use parkinds, only: rf=>kind_rf, im=>kind_im
use phys_fluid
use scm_io_default
#ifdef USE_NETCDF
use io, only: variable_metadata, write_2d_real_nc, write_1d_real_nc
#endif
use scm_sfx_data, only: taux, tauy, umst, hflx, qflx, cu
use scm_state_data
use pbl_turb_data
use pbl_turb_common
use pbl_dry_contrgradient
use pbl_solver
use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z, &
get_sigma_from_z_theta, get_theta_v, get_density_theta_vector
use diag_pbl
use pbl_grid
use sfx_data, only: meteoDataType, sfxDataType
use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
numericsType_most => numericsType
use config_parser
use config_utils
implicit none
!local functions
real, external:: get_ts_gabls2
!local varaibles
real time_begin, time_end, time_current
real dt, out_dt, output_dt
integer status, ierr
logical bstatus
character(len = 160) ::fname
character(len = 128) :: fname_nc
character(len = 128) ::fname_config = '../config-examples/config-gabls2-ex.txt'
real, allocatable :: ttemp(:,:), utemp(:,:),vtemp(:,:), ttime(:)
real, allocatable:: w_sub(:)
real(kind=rf):: ug, vg
real(kind=rf):: z0, zh, f_cor
real(kind=rf):: tsurf
real(kind=rf):: shir
type(fluidParamsDataType) fluid_params
type(pblgridDataType) grid
type(stateBLDataType):: bl, bl_old;
type(pblContrGradDataType):: cbl;
type(turbBLDataType):: turb
type(meteoDataType) :: meteo_cell
type(sfxDataType) :: sfx_cell
type(numericsType_most) :: numerics_sfx
!diagnostic variables
real hpbl
integer nkpbl
!io variables
real, dimension (5):: series_val
type (io_struct) :: series_f, scan_f
#ifdef USE_NETCDF
type(variable_metadata) :: meta_z, meta_t
type(variable_metadata) :: meta_temperature, meta_uwind, meta_vwind
#endif
integer k, nt, kmax
! init experiment params
numerics_sfx%maxiters_charnock = 10
ug = 3.00000
vg = -9.0
dt = 10.0
! call set_fluid_default(fluid_params)
! fluid_params%tref = 265.0
! fluid_params%pref = 1013.2500
! fluid_params%beta = fluid_params%g / fluid_params%tref
! fluid_params%kappa = 0.4
! fluid_params%p0= fluid_params%pref
!set output filenames
fname = 'test_gabls2.dat'
fname_nc = 'gabls2.nc'
series_f%fname = 'phys2.dat'
call set_file(series_f, series_f%fname)
scan_f%fname='time_scan2.dat'
call set_file(scan_f, scan_f%fname)
#ifdef USE_NETCDF
call set_meta_nc(meta_z, meta_t, meta_temperature, meta_uwind, meta_vwind)
#endif
call init_config(fname_config,status, ierr)
if (status == 0) then
write(*, *) ' FAILURE! > could not initialize config file: ', fname_config
ierr = 1 ! signal ERROR
stop
end if
call c_config_is_varname("time.begin"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("time.begin"//C_NULL_CHAR, time_begin, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: time.begin '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("time.end"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("time.end"//C_NULL_CHAR, time_end, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: time.end '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("time.dt"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("time.dt"//C_NULL_CHAR, dt, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: time.dt '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("time.out_dt"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("output.dt"//C_NULL_CHAR, output_dt, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: output.dt '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("Ug"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("Ug"//C_NULL_CHAR, ug, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: Ug '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("Vg"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("Vg"//C_NULL_CHAR, vg, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: Vg '
ierr = 1 ! signal ERROR
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
write(*, *) ' FAILURE! > could not set: Fluid params '
ierr = 1 ! signal ERROR
stop
end if
!z coordinate
call get_grid_params(grid, status, ierr)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: grid '
ierr = 1 ! signal ERROR
stop
end if
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)
bl_old%cor= f_cor
bl%cor = f_cor
bl%land_type = 0.0
bl%p0=fluid_params%p0
!setup surface
z0 = 0.03
zh = z0 / 10.0
! init blData states
call scm_data_allocate(bl, kmax)
call scm_data_allocate(bl_old, kmax)
call cbl_allocate(cbl, kmax)
call init_state(bl, ug, vg, fluid_params%tref)
call init_state(bl_old, ug, vg, fluid_params%tref)
call init_theta_gabls2(bl%theta,grid%z_cell, grid%kmax)
call scm_data_copy(bl_old,bl)
allocate(w_sub(0:kmax),source=0.0)
!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%kmax)
call get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
call scm_data_copy(bl_old,bl)
!Initialize turbulent closure
call turb_data_allocate(turb, kmax)
! getclosurefromconfig
call get_closure_params(turb, status, ierr)
!io setup
nt = floor((time_end - time_begin)/output_dt)
allocate(ttime(nt))
allocate(ttemp(kmax,nt),utemp(kmax,nt),vtemp(kmax,nt))
#ifdef USE_NETCDF
call write_1d_real_nc(grid%z_cell,fname_nc,meta_z)
#endif
out_dt = 0.0
nt=0
write(*,*)'nt begin', nt
call to_file_1d_2var('temperature2.dat', grid%z_cell, bl%theta, grid%kmax)
do while (time_current < time_end)
call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
tsurf = get_ts_gabls2(time_current)
meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%u(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
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
call get_density_theta_vector( bl, fluid_params, grid)
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) &
* 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
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)
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 add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
write(*,*) 'aftercoriolis,'
call get_wsub(w_sub, time_current, grid)
!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)
time_current = time_current + dt
if (time_current >=out_dt) then
nt = nt+1
ttime(nt) = time_current/3600.0
out_dt = out_dt + output_dt
write(*,*)'nt= ', nt
call theta2ta(bl_old%temp, bl_old%theta, &
fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call diag_pblh_inmcm(grid%z_cell,bl_old%theta,shir,0.0,umst,bl_old%kmax,hpbl)
!create output
series_val(1) = time_current/3600.0
series_val(2) = hflx
series_val(3) = turb%ustr
series_val(4) = hpbl
series_val(5) = 0.0
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)
#ifdef USE_NETCDF
do k = 1, kl
ttemp(k,nt) = thold(k)
utemp(k,nt) = uold(k)
vtemp(k,nt) = vold(k)
end do
#endif
!write(*,*) ccld * 10000.0
end if
end do
#ifdef USE_NETCDF
write(*,*)'nt= ', nt
call write_1d_real_nc(ttime,fname_nc,meta_t)
write(*,*)'here'
call write_2d_real_nc(ttemp,fname_nc,meta_temperature)
call write_2d_real_nc(utemp,fname_nc,meta_uwind)
call write_2d_real_nc(vtemp,fname_nc,meta_vwind)
#endif
call ta2theta(bl_old%theta, bl_old%temp, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call to_file_1d_3var(fname, grid%z_cell, bl_old%u, bl_old%v, bl_old%kmax)
!call to_file_1d_3var('coeffs.dat',grid%z_cell, bl_old%km, bl_old%kh, bl_old%kmax)
call to_file_1d_2var('temperature2.dat', grid%z_cell, bl%theta, bl_old%kmax)
call close_file(series_f)
call close_file(scan_f)
call deallocate_pbl_grid(grid)
call scm_data_deallocate(bl)
call scm_data_deallocate(bl_old)
call pbl_data_deallocate(turb)
call cbl_deallocate(cbl)
end program gabls2
subroutine init_state(bl, ug_,vg_,tref)
use parkinds, only: rf=>kind_rf, im=>kind_im
use scm_state_data, only:stateBLDataType
implicit none
real(kind=rf):: ug_,vg_,tref
type(stateBLDataType), intent(inout) :: bl
integer kmax
kmax= bl%kmax
write(*,*) 'blkmax=', kmax
bl%u(1:bl%kmax) = ug_
bl%v(1:bl%kmax) = vg_
bl%temp(1:bl%kmax) = tref
bl%theta(1:bl%kmax) = tref
bl%qv(1:bl%kmax) = 0.0
bl%lwp(1:bl%kmax) = 0.0
bl%cloud_frac(1:bl%kmax) = 0.0
bl%s_e(1:bl%kmax) = 0.0
bl%km(1:bl%kmax) = 0.0
bl%kh(1:bl%kmax) = 0.0
bl%vdcuv(1:bl%kmax) = 0.0
bl%vdctq(1:bl%kmax) = 0.0
end subroutine init_state
subroutine init_theta_gabls2(thold, z, kmax)
implicit none
integer,intent(in):: kmax
real, dimension(kmax):: thold, z
integer k
real T1, T2, Z1, Z2
!write(*,*) 'theta_init'
do k = kmax, 1, -1
if (z(k) <= 200.0 ) then
z1 = 0.0
z2 = 200.0
T1 = 288.0
T2 = 286.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 200.0 .and. z(k)<= 850.0 ) then
z1 = 200.0
z2 = 850.0
T1 = 286.0
T2 = 286.0
thold(k) = t2
elseif (z(k) > 850.0 .and. z(k)<= 900.0) then
z1 = 850.0
z2 = 900.0
T1 = 286.0
T2 = 288.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 900.0 .and. z(k)<= 1000.0) then
z1 = 900.0
z2 = 1000.0
T1 = 288.0
T2 = 292.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 1000.0 .and. z(k)<= 2000.0) then
z1 = 1000.0
z2 = 2000.0
T1 = 292.0
T2 = 300.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 2000.0 .and. z(k)<= 3500.0) then
z1 = 2000.0
z2 = 3500.0
T1 = 300.0
T2 = 310.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 3500.0 .and. z(k)<= 4000.0) then
z1 = 3500.0
z2 = 4000.0
T1 = 310.0
T2 = 312.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 4000.0) then
thold(k) = 312.0
end if
!write(*,*) z(k), thold(k)
end do
end subroutine init_theta_gabls2
real function get_ts_gabls2(curtime)
implicit none
real curtime
real hrs, ts
hrs = curtime / 3600.0 + 16.0
ts = 4.4
if (hrs <= 17.4) then
ts = -10.0 - 25.0 * cos(0.22 * hrs + 0.2)
end if
if ((hrs > 17.4) .and. (hrs <= 30.0)) then
ts = -0.54 * hrs + 15.2
end if
if ((hrs > 30.0) .and. (hrs <= 41.9)) then
ts = -7.0 - 25.0 * cos(0.21 * hrs + 1.8)
end if
if ((hrs > 41.9) .and. (hrs <= 53.3)) then
ts = -0.37 * hrs + 18.0;
end if
if ((hrs > 53.3) .and. (hrs <= 65.6)) then
ts = -4.0 - 25.0 * cos(0.22 * hrs + 2.5)
end if
ts = ts + 273.15
get_ts_gabls2 = ts
end function get_ts_gabls2
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(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)
end do
end subroutine add_coriolis
subroutine get_wsub(w, curtime, grid)
use pbl_grid, only: pblgridDataType
implicit none
type(pblgridDataType), intent(in):: grid
real, intent(in):: curtime
real, dimension(*), intent(out):: w
integer k
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)
else
w(k) = -0.005
end if
end do
else
w(1:grid%kmax) = 0.0
end if
end subroutine get_wsub
subroutine put_tscan( time, z, bl, nl, f)
use parkinds, only : rf=>kind_rf, im => kind_im
use scm_io_default
use scm_state_data, only:stateBLDataType
implicit none
type(stateBLDataType), intent(in):: bl
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) = bl%theta(k)
stamp(4,k) = bl%u(k)
stamp(5,k) = bl%v(k)
end do
! call to write timestamp
call write_timescan(stamp,nl, 5, f)
end subroutine put_tscan
subroutine get_surface_from_config(model, surface, z0,zh )
use iso_c_binding, only: C_NULL_CHAR
use config_parser
use sfx_common, only: char_array2str
use sfx_config
use sfx_surface, only: get_surface_id
use parkinds, only: rf=>kind_rf
implicit none
character, allocatable :: config_field(:)
integer,intent(out) :: model, surface
real(kind=rf), intent(out) :: z0,zh
integer status, ierr
call c_config_is_varname("model.id"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_string("model.id"//C_NULL_CHAR, config_field, status)
if (status == 0) then
ierr = 1 ! signal ERROR
stop
end if
model = get_model_id(char_array2str(config_field))
if (model == -1) then
write(*, *) ' FAILURE! > unknown model [key]: ', trim(char_array2str(config_field))
ierr = 1 ! signal ERROR
return
end if
end if
call c_config_is_varname("surface.type"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_string("surface.type"//C_NULL_CHAR, config_field, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
surface = get_surface_id(char_array2str(config_field))
if (surface == -1) then
write(*, *) ' FAILURE! > unknown surface [key]: ', trim(char_array2str(config_field))
ierr = 1 ! signal ERROR
return
end if
endif
end subroutine get_surface_from_config
#ifdef USE_NETCDF
subroutine set_meta_nc(z_meta, t_meta, theta_meta, uwind_meta, vwind_meta)
use io, only: variable_metadata
type(variable_metadata), intent(inout):: z_meta, t_meta, theta_meta, uwind_meta, vwind_meta
z_meta = variable_metadata( &
name = 'Z', &
dim_names = [character(len=32) :: 'Z', '', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Height', &
'm', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
t_meta = variable_metadata( &
name = 'Time', &
dim_names = [character(len=32) :: 'Time', '', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Time', &
'hours since 00-00-00', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
theta_meta = variable_metadata( &
name = 'th', &
dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Potential Temperature', &
'K', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
uwind_meta = variable_metadata( &
name = 'ua', &
dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Longitude wind component', &
'm/s', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
vwind_meta = variable_metadata( &
name = 'va', &
dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
char_name = [character(len=32) :: &
'long_name', &
'units', &
'', '', '', '', '', '', '', ''], &
char_value = [character(len=32) :: &
'Latitude wind component', &
'm/s', &
'', '', '', '', '', '', '', ''], &
real_name = [character(len=32) :: &
'missing_value', &
'scale_factor', &
'', '', '', '', '', '', '', ''], &
real_value = [&
-9999.9, &
1.0, &
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
end subroutine set_meta_nc
#endif
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment