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

reworked some functions

parent a914e6d8
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 640 additions and 91 deletions
CMakeLists.txt
+ 7
− 2
View file @ ed08e8fc
......@@ -36,6 +36,9 @@ set(lib_files
src/scm_state_data.f90
src/pbl_turb_data.f90
src/pbl_solver.f90
src/fo_stability_functions.f90
src/pbl_fo_turb.f90
src/pbl_turb_common.f90
src/diag_pbl.f90)
add_library(abl_lib ${lib_files})
......@@ -117,4 +120,6 @@ if (WITH_TESTS)
PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -fbacktrace -ffpe-trap=zero,overflow,underflow >)
endif()
endif()
endif()
config-examples/config-gabls1-ex.txt
+ 25
− 9
View file @ ed08e8fc
Ug = 8.0;
Vg = 0.0;
hbl_0 = 100.0;
Tgrad = 0.01;
cooling_rate = -0.25 * 3600.0;
time{
begin = 0.0;
end = 9.0 * 3600.0;
end = 20.0; #9.0 * 3600.0;
dt = 20.0;
out_dt = 360.0;
}
......@@ -16,8 +20,15 @@ fluid {
kappa = 0.4; # von Karman constant
}
surface_model {
grid {
type = "uniform";
nz = 8;
h_bot = 0.0;
h_top = 400.0;
}
surface_model {
id = "esm"; # optional: "esm" (default), "sheba", "most", "log"
type = "land";
z0_m = 0.01; # aerodynamic roughness [m]
......@@ -27,3 +38,8 @@ surface_model {
output{
dt = 0.2 *3600.0;
}
closure {
type = "FOM";
name = "Louis"; # "Esau"; # "EFB";
}
\ No newline at end of file
src/config-utils.f90
+ 64
− 0
View file @ ed08e8fc
......@@ -8,6 +8,7 @@ module config_utils
implicit none
integer, public, save:: is_config_initialized = 0
public init_config, get_fluid_params, get_grid_params
public get_closure_params
!public get_geo_forcing, get_heat_forcing
......@@ -86,6 +87,7 @@ module config_utils
grid_inmcm21_tag, grid_inmcm73_tag, &
grid_streached_tag, grid_uniform_tag, &
set_pbl_grid_uniform
use sfx_common, only: char_array2str
implicit none
......@@ -105,15 +107,18 @@ module config_utils
if ((status /= 0)) then
call c_config_get_string("grid.type"//C_NULL_CHAR, config_field, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: grid.type '
ierr = 1 ! signal ERROR
return
end if
end if
tag = char_array2str(config_field)
if (trim(tag) == trim(grid_uniform_tag)) then
call c_config_is_varname("grid.nz"//C_NULL_CHAR, status)
if ((status /= 0)) then
call c_config_get_int("grid.nz"//C_NULL_CHAR, nz, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: grid.nz '
ierr = 1 ! signal ERROR
return
end if
......@@ -123,6 +128,7 @@ module config_utils
if ((status /= 0)) then
call c_config_get_float("grid.h_bot"//C_NULL_CHAR, h_bot, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: grid.h_bot '
ierr = 1 ! signal ERROR
return
end if
......@@ -131,6 +137,7 @@ module config_utils
if ((status /= 0)) then
call c_config_get_float("grid.h_top"//C_NULL_CHAR, h_top, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: grid.h_top '
ierr = 1 ! signal ERROR
return
end if
......@@ -144,6 +151,63 @@ module config_utils
end subroutine get_grid_params
subroutine get_closure_params(turb, status, ierr)
use pbl_turb_data, only: turbBLDataType, &
pbl_fom_tag, pbl_kl_tag, &
pbl_fom_louis_tag, pbl_fom_esau_tag, &
pbl_fom_efb_tag
use sfx_common, only: char_array2str
implicit none
type(turbBLDataType), intent(inout):: turb
integer,intent(out):: status, ierr
character(len=50) :: tag, tag_stabf
character, allocatable :: config_field(:)
ierr = 0
if( is_config_initialized /= 0 ) then
! closure type
call c_config_is_varname("closure.type"//C_NULL_CHAR, status)
if ((status /= 0)) then
call c_config_get_string("closure.type"//C_NULL_CHAR, config_field, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: grid.type '
ierr = 1 ! signal ERROR
return
end if
end if
tag = char_array2str(config_field)
if (trim(tag) == trim(pbl_fom_tag)) then
turb%cl_type = 1
call c_config_is_varname("closure.name"//C_NULL_CHAR, status)
if ((status /= 0)) then
call c_config_get_string("closure.name"//C_NULL_CHAR, config_field, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: closure.name '
ierr = 1 ! signal ERROR
return
else
tag_stabf = char_array2str(config_field)
if (trim(tag_stabf) == trim(pbl_fom_louis_tag)) then
turb%sf_type = 1
elseif (trim(tag_stabf) == trim(pbl_fom_esau_tag)) then
turb%sf_type = 2
elseif (trim(tag_stabf) == trim(pbl_fom_efb_tag)) then
turb%sf_type = 3
else
status = 1
ierr = 1
write(*, *) ' FAILURE! > could not find eligable: closure.name '
write(*,*) tag_stabf, pbl_fom_louis_tag
return
end if
end if
end if
else
turb%cl_type = 2 ! KL - type
endif
end if
end subroutine get_closure_params
!> @brief character array to string conversion
function char_array2str(char_array) result(str)
......
src/fo_stability_functions.f90 0 → 100644
+ 91
− 0
View file @ ed08e8fc
! Created by Andrey Debolskiy on 24.01.2025.
module fo_stability_functions
!< Constants to for stability functions
!< louis stability functions
real, parameter:: yb = 5.e0
real, parameter:: yc = 5.e0
real, parameter:: yd = 5.e0
real, parameter:: ye = yb*yc*sqrt(1.e0/3.e0)
real, parameter::ricr = 2.e0/(3.e0*yd)
real, parameter:: vlinfm = 120.0e0
real, parameter:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)
! security parameters
! du2min is a minimal value to calculate wind shear
! 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
contains
!>
subroutine louis_stab_functions(fnriuv,fnritq,rinum,vdlm, &
& vdlh,zkh, rolim)
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
if(rinum.ge.0.e0) then
scf = sqrt(1.e0 + yd*rinum)
fnriuv = 1.e0/(1.e0 + 2.e0*yb * rinum / scf)
fnritq = 1.e0/(1.e0 + 3.e0*yb * rinum * scf)
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(rolim.gt.0.99.and.rolim.lt.1.5) 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 louis_stab_functions
!>
!>esau byrkjedal 2007 + viterbo 1999 unstable stratification
subroutine esau_vit_stab_functions( fnriuv,fnritq,rinum, &
& vdlm,zkh, dz)
implicit none
real fnriuv,fnritq,rinum, vdlm, zkh, dz
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
gh = ((vdlm * vdlm)/(dz * sqrt(abs(dz * zkh)))) &
& * ((1 + dz / zkh)**(1/3) - 1)**(3/2)
fnriuv = 1.0 - (2.0 * bb * rinum) &
& / (1.0 + 3.0 * bb * cc * gh * sqrt(abs(rinum)))
fnritq = 1.0 - (2.0 * bb * rinum) &
& / (1.0 + 3.0 * bb * cc * gh * sqrt(abs(rinum)))
end if
end subroutine esau_vit_stab_functions
end module fo_stability_functions
\ No newline at end of file
src/pbl_fo_turb.f90
+ 105
− 16
View file @ ed08e8fc
......@@ -2,25 +2,93 @@
module pbl_fo_turb
implicit none
!constants to for stability functions
real, parameter:: yb = 5.e0
real, parameter:: yc = 5.e0
real, parameter:: yd = 5.e0
real:: ye = yb*yc*sqrt(1.e0/3.e0)
real:: ricr = 2.e0/(3.e0*yd)
real, parameter:: vlinfm = 160.0e0
real:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)
real, parameter:: aka = 0.4e0
!real, parameter:: r = 287.05
!real, parameter:: cp = 1.005
! security parameters
! du2min is a minimal value to calculate wind shear
! 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
public get_fo_coeffs
public get_fo_coeffs_louis
public get_turb_length
public get_fo_diff
public get_rig
! public :: default_stab_functions
! public :: esau_stab_functions
! private:: get_ph_efb
! private:: q_root
contains
subroutine get_fo_diff(turb, bl, fluid, grid)
! subroutine get_fo_diff(turb, bl, fluid, grid)
! use pbl_turb_data, only: turbBLDataType
! use scm_state_data, only: stateBLDataType
! use pbl_grid, only: pblgridDataType
! use phys_fluid, only: fluidParamsDataType
! end subroutine get_fo_diff
subroutine get_fo_coeffs(turb, bl, fluid, grid)
use pbl_turb_data, only: turbBLDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
use phys_fluid, only: fluidParamsDataType
end subroutine get_fo_diff
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(inout):: bl
type(pblgridDataType), intent(in)::grid
type(turbBLDataType), intent(inout):: turb
subroutine get_fo_coeffs(turb, bl, fluid, grid)
if (turb%sf_type == 1) then !> Louis stab functions
call get_fo_coeffs_louis(turb, bl, fluid, grid)
elseif (turb%sf_type == 2) then !> Esau stab functions
! call get_fo_coeffs_esau(turb, bl, fluid, grid)
elseif (turb%sf_type == 3) then !> EFB stab functions
! call get_fo_coeffs_efb(turb, bl, fluid, grid)
else
write(*,*) 'FAILURE: stab functions type in FOM not set!'
end if
end subroutine get_fo_coeffs
subroutine get_fo_coeffs_louis(turb, bl, fluid, grid)
use pbl_turb_data, only: turbBLDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
use phys_fluid, only: fluidParamsDataType
end subroutine get_fo_coeffs
use fo_stability_functions
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(inout):: bl
type(pblgridDataType), intent(in)::grid
type(turbBLDataType), intent(inout):: turb
real fnriuv, fnritq,rho
real ocean_flag
integer k, kmax
ocean_flag = real(bl%land_type)
do k = kmax, 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)
end do
! for compliance with old version
bl%vdcuv(grid%kmax) = 0.0
bl%vdcuv(grid%kmax) = 0.0
do k = kmax-1,1, -1
bl%vdcuv(k) = 0.5 * (bl%rho(k) + bl%rho(k+1)) * turb%km(k)
bl%vdctq(k) = 0.5 * (bl%rho(k) + bl%rho(k+1)) * turb%kh(k)
end do
end subroutine get_fo_coeffs_louis
subroutine get_turb_length(turb, bl, fluid, grid, hbl_option)
use pbl_turb_data, only: turbBLDataType
......@@ -30,19 +98,40 @@ module pbl_fo_turb
use diag_pbl
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(in):: bl
type(stateBLDataType), intent(inout):: bl
type(pblgridDataType), intent(in)::grid
type(turbBLDataTyoe), intent(inout):: turb
type(turbBLDataType), intent(inout):: turb
integer hbl_option
real l_neutral, lm, lh
call get_rig(turb, bl, fluid, grid)
integer k, 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)
else
write(*,*) 'wrong hbl diagnostics option'
return
end if
do k = kmax,1, -1
l_neutral = fluid%kappa * grid%z_cell(k)
if ( k <= bl%kpbl) then
lm = 0.0
lh = 0.0
else
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
lm = 1.0
lh = 1.0
endif
endif
turb%l_turb_m = l_neutral * lm
turb%l_turb_h = l_neutral * lh
end do
end subroutine get_turb_length
subroutine get_rig(turb, bl, fluid, grid)
use pbl_turb_data, only: turbBLDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
use phys_fluid, only: fluidParamsDataType
end subroutine get_rig
end module pbl_fo_turb
\ No newline at end of file
src/pbl_grid.f90
+ 1
− 0
View file @ ed08e8fc
......@@ -95,6 +95,7 @@ module pbl_grid
call allocate_pbl_grid(grid, nk)
dz = (h_top - h_bottom) / real(nk, kind = kind(dz))
grid%z_edge(nk) = h_bottom
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))
......
src/pbl_turb_common.f90 0 → 100644
+ 112
− 0
View file @ ed08e8fc
! Created by Andrey Debolskiy on 11.02.2025.
!> Collection of utilities for pbl_turb_data
module pbl_turb_common
use pbl_turb_data
implicit none
public get_s2, get_n2, get_rig
public get_coeffs_general
contains
subroutine get_n2(turb, bl, fluid, grid)
use pbl_turb_data, only: turbBLDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
use phys_fluid, only: fluidParamsDataType
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(in):: bl
type(pblgridDataType), intent(in)::grid
type(turbBLDataType), intent(inout):: turb
real dtvir, buoyp
integer k
do k = 1,grid%kmax-1
dtvir = bl%theta_v(k) - bl%theta_v(k+1)
buoyp = fluid%g / ( (bl%theta_v(k+1) + bl%theta_v(k))/2.0e0 )
turb%n2 = buoyp*(grid%z_cell(k) - grid%z_cell(k+1))*dtvir
end do
end subroutine get_n2
subroutine get_s2(turb, bl, fluid, grid, du2min)
use pbl_turb_data, only: turbBLDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
use phys_fluid, only: fluidParamsDataType
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(in):: bl
type(pblgridDataType), intent(in)::grid
type(turbBLDataType), intent(inout):: turb
real, intent(in) :: du2min
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)) &
* (grid%z_cell(k) - grid%z_cell(k+1)))
enddo
end subroutine get_s2
subroutine get_rig(turb, bl, fluid, grid, du2min, ricr)
use pbl_turb_data, only: turbBLDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
use phys_fluid, only: fluidParamsDataType
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(in):: bl
type(pblgridDataType), intent(in)::grid
real, intent(in):: du2min, ricr
type(turbBLDataType), intent(inout):: turb
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)
buoyp = fluid%g / ( (bl%theta_v(k+1) + bl%theta_v(k))/2.0e0 )
rinum = buoyp*(grid%z_cell(k) - grid%z_cell(k+1))*dtvir/du2
rinum = amin1(ricr,rinum)
turb%rig(k) = rinum
end do
end subroutine get_rig
subroutine get_coeffs_general(turb, bl, fluid, hbl_option, grid)
use pbl_turb_data, only: turbBLDataType
use scm_state_data, only: stateBLDataType
use pbl_grid, only: pblgridDataType
use phys_fluid, only: fluidParamsDataType
use pbl_fo_turb, only: get_fo_coeffs, get_turb_length, du2min, ricr
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(inout):: bl
type(pblgridDataType), intent(in)::grid
type(turbBLDataType), intent(inout):: turb
integer, intent(in):: hbl_option
write(*,*) 'here1'
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_turb_length(turb, bl, fluid, grid, hbl_option)
if (turb%cl_type == 1) then ! FOM closures
write(*,*) 'here4'
call get_fo_coeffs(turb, bl, fluid, grid)
end if
if (turb%cl_type == 2) then ! KL closures
!call get_kl_coeffs(turb, bl, fluid, grid)
end if
end subroutine get_coeffs_general
end module pbl_turb_common
\ No newline at end of file
src/pbl_turb_data.f90
+ 40
− 12
View file @ ed08e8fc
! Created by Andrey Debolskiy on 26.11.2024.
module pbl_turb_data
!<@brief turbulent fields
type turbBLDataType
real, allocatable:: uw(:), vw(:), thetaw(:), thetavw(:), qvw(:)
real, allocatable:: s2(:),n2(:)
real, allocatable :: qcw(:), cloud_frac_w(:)
real, allocatable:: tke(:), eps(:), l_turb(:)
real, allocatable:: rig(:), prt(:), rif(:)
real, allocatable:: km(:), kh(:)
real, allocatable:: uw(:) !< u-component of kinematic vertical momentum flux, [m2/s2]
real, allocatable:: vw(:) !< v-component of kinematic vertical momentum flux, [m2/s2]
real, allocatable:: thetaw(:) !< kinematic vertical heat flux flux, [K m/s2]
real, allocatable::thetavw(:) !< kinematic vertical virtual temperature flux, [ K m/s2]
real, allocatable:: qvw(:) !< kinematic vertical moisture flux, [ bar m/s]
real, allocatable:: s2(:) !< square of shear frequancy [s-2]
real, allocatable:: n2(:) !< Brunt-Vaisala frequancy [s-2]
real, allocatable :: qcw(:) !< turbulent vertical liquid water flux [kg/kg m /s2]
real, allocatable:: cloud_frac_w(:) !< turbulent vertical cloud fraction flux [ m /s2]
real, allocatable:: tke(:) !< turbulent kinetic energy [m2/s2]
real, allocatable:: eps(:) !< dissipation rate of turbulent kinetic energy [m2/s3]
real, allocatable:: l_turb_m(:) !< turbulent length scale for momentum [m]
real, allocatable:: l_turb_h(:) !< turbulent length scale for heat [m]
real, allocatable:: rig(:) !< gradient Richardson number
real, allocatable::prt(:) !< turbulent Prandtl number
real, allocatable::rif(:) !< flux Richardson number
real, allocatable:: km(:) !< kinematic turbulent viscosity coefficient [ m2 / s]
real, allocatable:: kh(:) !< kinematic turbulent diffusion coefficient [ m2 / s]
real:: ustr !< surface dynamic velocity scale
real:: tstar !< surface dynamic temperature scale
integer :: kmax
integer :: cl_type !< closure type - "FOM"/"KL"
integer :: sf_type !< stability functions type - "Louis"/"Esau"/"EFB"
end type turbBLDataType
public turb_data_allocate, turb_data_deallocate
!closure types
character(len = 16), parameter :: pbl_fom_tag = 'FOM'
character(len = 16), parameter :: pbl_kl_tag = 'KL'
!FOM stab functions
character(len = 16), parameter :: pbl_fom_louis_tag = 'Louis'
character(len = 16), parameter :: pbl_fom_esau_tag= 'Esau'
character(len = 16), parameter :: pbl_fom_efb_tag = 'EFB'
public turb_data_allocate
public turb_data_deallocate
public update_fields
contains
subroutine turb_data_allocate(turb_data, kmax)
implicit none
......@@ -29,7 +56,8 @@ module pbl_turb_data
allocate(turb_data%n2(kmax))
allocate(turb_data%tke(kmax))
allocate(turb_data%eps(kmax))
allocate(turb_data%l_turb(kmax))
allocate(turb_data%l_turb_h(kmax))
allocate(turb_data%l_turb_m(kmax))
allocate(turb_data%rig(kmax))
allocate(turb_data%rif(kmax))
allocate(turb_data%prt(kmax))
......@@ -37,7 +65,7 @@ module pbl_turb_data
allocate(turb_data%kh(kmax))
end subroutine turb_data_allocate
subroutine scm_data_deallocate(turb_data)
subroutine pbl_data_deallocate(turb_data)
implicit none
type(turbBLDataType), intent(inout):: turb_data
deallocate(turb_data%uw)
......@@ -53,10 +81,10 @@ module pbl_turb_data
deallocate(turb_data%kh)
deallocate(turb_data%tke)
deallocate(turb_data%eps)
deallocate(turb_data%l_turb)
deallocate(turb_data%l_turb_m)
deallocate(turb_data%l_turb_h)
deallocate(turb_data%rig)
deallocate(turb_data%rif)
deallocate(turb_data%prt)
end subroutine scm_data_deallocate
end subroutine pbl_data_deallocate
end module pbl_turb_data
\ No newline at end of file
src/scm_state_data.f90
+ 31
− 18
View file @ ed08e8fc
! Created by Andrey Debolskiy on 20.11.2024.
module scm_state_data
!< @brief abl state def.
implicit none
type stateBLDataType
real, allocatable:: u(:), v(:), temp(:), theta(:) ,qv(:)
real, allocatable:: rho(:)
real, allocatable :: lwp(:), cloud_frac(:), s_e(:)
real, allocatable:: km(:), kh(:)
real, allocatable:: vdctq(:)
real, allocatable:: vdcuv(:)
real :: HPBLA_diag
real :: p0
integer :: ktvdm=10
integer :: ktvd
integer :: kpbl
integer :: ktop
integer :: kmax
real, allocatable:: u(:) !< u-componet velocity, [m/s]
real, allocatable:: v(:) !< v-componet velocity, [m/s]
real, allocatable:: temp(:) !< absolute temperature, [K]
real, allocatable:: theta(:) !< potential temperature, [K]
real, allocatable:: qv(:) ! specific humidity,[kg/m3]
real, allocatable:: theta_v(:) !< virtual potential temperature, [K]
real, allocatable:: rho(:) !< air density, [K]
real, allocatable:: lwp(:) !< liquid water path [kg/kg]
real, allocatable:: cloud_frac(:) !< cloud fraction [0-1]
real, allocatable:: s_e(:) !< Static energy [ m2/ s2]
real, allocatable:: km(:) !< kinematic viscosity koefficient [m2 / s2]
real, allocatable:: kh(:) !< kinematic diffusivity koefficient [ K m / s2]
real, allocatable:: vdcuv(:) !< viscosity koefficient [kg /m s2]
real, allocatable:: vdctq(:) !< diffusivity koefficient [ kg K / m2 s2]
real :: p0 !< surface pressure [bar]
real :: cor !< in column corioulis parameter [s-1]
integer :: kmax !< number of vertical cells
real :: hpbla_diag !< diagnostic boundary layer height [m]
integer :: ktvdm !< top layer index for diffusion
integer :: kpbl !< boundary layer height index
integer :: land_type !< ocean/land flag
end type stateBLDataType
public scm_data_allocate
public scm_data_deallocate
......@@ -30,14 +40,15 @@ module scm_state_data
allocate(bl_data%v(kmax))
allocate(bl_data%temp(kmax))
allocate(bl_data%theta(kmax))
allocate(bl_data%theta_v(kmax))
allocate(bl_data%qv(kmax))
allocate(bl_data%lwp(kmax))
allocate(bl_data%cloud_frac(kmax))
allocate(bl_data%s_e(kmax))
allocate(bl_data%km(kmax))
allocate(bl_data%kh(kmax))
allocate(bl_data%vdcuv(kmax))
allocate(bl_data%vdctq(kmax))
allocate(bl_data%vdcuv(0:kmax))
allocate(bl_data%vdctq(0:kmax))
end subroutine scm_data_allocate
subroutine scm_data_deallocate(bl_data)
......@@ -47,6 +58,7 @@ module scm_state_data
deallocate(bl_data%v)
deallocate(bl_data%temp)
deallocate(bl_data%theta)
deallocate(bl_data%theta_v)
deallocate(bl_data%qv)
deallocate(bl_data%lwp)
deallocate(bl_data%cloud_frac)
......@@ -60,14 +72,15 @@ module scm_state_data
subroutine scm_data_copy(bl, bl_old)
implicit none
type(stateBLDataType),intent(in):: bl_old
type(stateBLDataType),intent(out):: bl
type(stateBLDataType),intent(inout):: bl
if (bl%kmax /= bl_old%kmax) then
write(*,*) 'error in copy BLData: size is not compatible'
write(*,*) 'error in copy BLData: size is not compatible', bl%kmax, bl_old%kmax
else
bl%u = bl_old%u
bl%v = bl_old%v
bl%temp = bl_old%temp
bl%theta = bl_old%theta
bl%theta_v = bl_old%theta_v
bl%qv = bl_old%qv
bl%lwp = bl_old%lwp
bl%cloud_frac= bl_old%cloud_frac
......@@ -80,7 +93,7 @@ module scm_state_data
bl%p0 = bl_old%p0
bl%ktvdm = bl_old%ktvdm
bl%kpbl = bl_old%kpbl
bl%ktop = bl_old%ktop
bl%cor = bl_old%cor
end if
end subroutine scm_data_copy
......
src/state_utils.f90
+ 64
− 3
View file @ ed08e8fc
......@@ -9,10 +9,12 @@ module state_utils
public get_exner
public get_sigma_from_z
public get_z_from_sigma
public get_z_from_sigma_theta
public get_density_theta
public get_density
public dqsdt_sc, qmax_sc, qmax
public DQSDT2, DQSDT, CLDT
public get_theta_v
contains
subroutine get_geopotential_moist(F,t, qv, sig, fs, rd, eps, kmax)
implicit none
......@@ -186,6 +188,22 @@ module state_utils
get_density_theta = (p / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
end function get_density_theta
subroutine get_theta_v( theta_v, theta, qv, fluid, kmax)
use phys_fluid, only: fluidParamsDataType
implicit none
integer,intent(in):: kmax
type(fluidParamsDataType), intent(in):: fluid
real,dimension(kmax), intent(in):: theta
real,dimension(kmax), intent(in):: qv
real,dimension(kmax), intent(out):: theta_v
integer k
do k = kmax,1,-1
theta_v(k) = theta(k) * (1.0 + fluid%R_d / fluid%R_v * qv(k))
end do
end subroutine get_theta_v
subroutine get_exner( exnr, p, fluid, kmax)
use phys_fluid, only: fluidParamsDataType
implicit none
......@@ -233,10 +251,12 @@ module state_utils
kmax = grid%kmax
grid%sigma_edge(kmax) = 1.0
rho = get_density(ps, bl%temp(kmax), bl%qv(kmax), fluid)
dp = - rho * fluid%g * grid%dzc(kmax)
grid%sigma(kmax) = (ps + dp/2) / ps
grid%sigma_edge(kmax) = (ps + dp) / ps
grid%sigma(kmax) = 1.0 + (dp/2) / ps
grid%sigma_edge(kmax) = 1.0
p = ps + dp
do k = kmax-1,1,-1
rho = get_density(p, bl%temp(k), bl%qv(k), fluid)
dp = - rho * fluid%g * grid%dzc(k)
......@@ -244,7 +264,7 @@ module state_utils
rho = get_density(p, 0.5*(bl%temp(k+1) + bl%temp(k)), &
0.5*(bl%qv(k+1) + bl%qv(k)), fluid)
grid%sigma(k) = grid%sigma(k+1) &
- rho * fluid%g * grid%dze(k)
- rho * fluid%g * grid%dze(k) / ps
p = grid%sigma_edge(k) * ps
end do
grid%dsigmac(2:kmax) = grid%sigma(2:kmax) - grid%sigma(1:kmax-1)
......@@ -252,6 +272,47 @@ module state_utils
grid%dsigmae(1:kmax) = grid%sigma_edge(1:kmax) - grid%sigma_edge(1:kmax-1)
end subroutine get_sigma_from_z
subroutine get_sigma_from_z_theta( grid, bl, ps, fluid)
use pbl_grid, only : pblgridDataType
use phys_fluid, only: fluidParamsDataType
use scm_state_data, only: stateBLDataType
implicit none
type(fluidParamsDataType), intent(in):: fluid
type(stateBLDataType), intent(in):: bl
real, intent(in):: ps
type(pblgridDataType), intent(inout)::grid
real dp, rho, p, temp, appa
integer k, kmax
kmax = grid%kmax
grid%sigma_edge(kmax) = 1.0
!convert theta -> temp
appa = fluid%R_d / fluid%cp
temp = bl%theta(kmax) * (100000.0 / ps )**appa
rho = get_density(ps, temp, bl%qv(kmax), fluid)
dp = - rho * fluid%g * grid%dzc(kmax)
grid%sigma(kmax) = 1.0 + (dp/2) / ps
grid%sigma_edge(kmax) = 1.0
p = ps + dp
do k = kmax-1,1,-1
temp = bl%theta(k) * (100000.0 / p )**appa
rho = get_density(p, bl%temp(k), bl%qv(k), fluid)
dp = - rho * fluid%g * grid%dzc(k)
grid%sigma_edge(k) = grid%sigma_edge(k+1) + dp / ps
rho = get_density(p, 0.5*(bl%temp(k+1) + bl%temp(k)), &
0.5*(bl%qv(k+1) + bl%qv(k)), fluid)
grid%sigma(k) = grid%sigma(k+1) &
- rho * fluid%g * grid%dze(k) / ps
p = grid%sigma_edge(k) * ps
end do
grid%dsigmac(2:kmax) = grid%sigma(2:kmax) - grid%sigma(1:kmax-1)
grid%sigma_edge(0) = 0.0
grid%dsigmae(1:kmax) = grid%sigma_edge(1:kmax) - grid%sigma_edge(1:kmax-1)
end subroutine get_sigma_from_z_theta
subroutine get_z_from_sigma( grid, bl, zs, fluid, kmax)
use pbl_grid, only : pblgridDataType
use phys_fluid, only: fluidParamsDataType
......
src/tests/gabls1.f90
+ 100
− 31
View file @ ed08e8fc
......@@ -10,8 +10,10 @@ program gabls1
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_solver
use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z
use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z, &
get_sigma_from_z_theta, get_theta_v
use diag_pbl
use pbl_grid
use sfx_data, only: meteoDataType, sfxDataType
......@@ -35,6 +37,9 @@ program gabls1
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):: cooling_rate
type(fluidParamsDataType) fluid_params
type(pblgridDataType) grid
......@@ -89,8 +94,9 @@ program gabls1
call init_config(fname_config,status, ierr)
if (status == 0) then
write(*, *) ' FAILURE! > could not initialize config file: ', fname_config
ierr = 1 ! signal ERROR
return
stop
end if
call c_config_is_varname("time.begin"//C_NULL_CHAR, status)
......@@ -98,8 +104,9 @@ program gabls1
!< 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
return
stop
end if
end if
call c_config_is_varname("time.end"//C_NULL_CHAR, status)
......@@ -107,8 +114,9 @@ program gabls1
!< 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
return
stop
end if
end if
call c_config_is_varname("time.dt"//C_NULL_CHAR, status)
......@@ -116,18 +124,20 @@ program gabls1
!< 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
return
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("time.out_dt"//C_NULL_CHAR, output_dt, status)
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
return
stop
end if
end if
......@@ -136,8 +146,9 @@ program gabls1
!< 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
return
stop
end if
end if
call c_config_is_varname("Vg"//C_NULL_CHAR, status)
......@@ -145,33 +156,78 @@ program gabls1
!< 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
return
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
return
stop
end if
call c_config_is_varname("cooling_rate"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("cooling_rate"//C_NULL_CHAR, cooling_rate, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: cooling_rate '
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
return
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
write(*,*) 'F coriolis: = ', f_cor
!fs(1, 1) = 0.0
!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)
......@@ -179,22 +235,27 @@ program gabls1
!TODO why zet is not initialized in beirit properly?????
!CALL GEO(, FS(1, 1), grid%z_cell, grid%kmax) !change for clima
DO K = 1, kmax
grid%z_cell(k) = grid%z_cell(k) / fluid_params%g
DO K = 1, kmax
if (grid%z_cell(k) > 100.0) then
bl_old%theta(k) = fluid_params%tref + 0.01* (grid%z_cell(k) - 100.0)
bl_old%theta(k) = fluid_params%tref + tgrad * (grid%z_cell(k) - 100.0)
else
bl_old%theta(k) = fluid_params%tref
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)
bl_old%theta_v(1:kmax) = bl_old%theta(1:kmax) * (1.0+ fluid_params%eps_v * bl_old%qv(1:kmax))
call scm_data_copy(bl,bl_old)
!Initialize turbulent closure
call turb_data_allocate(turb, kmax)
! getclosurefromconfig
call get_closure_params(turb, status, ierr)
!finish updating grid
call get_sigma_from_z( grid, bl, fluid_params%p0 * 100.0, fluid_params)
!io setup
......@@ -208,7 +269,8 @@ program gabls1
nt=0
write(*,*)'nt=0', nt
do while (time_current < time_end)
tsurf = fluid_params%tref - 0.25 * (time_current - time_begin) / 3600.0;
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)
meteo_cell%dT = -tsurf + bl_old%theta(kmax)
meteo_cell%Tsemi = 0.5 * (tsurf + bl_old%theta(kmax))
......@@ -216,16 +278,20 @@ program gabls1
meteo_cell%h = grid%z_cell(kmax)
meteo_cell%z0_m = z0
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
cu = 1000.0 * sfx_cell%Cm**2
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
umst = sfx_cell%Cm * meteo_cell%U
turb%ustr = sfx_cell%Cm * meteo_cell%U
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)
!call add_coriolis(bl, ug, vg, f_cor, dt, grid)
call scm_data_copy(bl,bl_old)
call get_fo_diff(turb, bl, grid)
call solv_diffusion(bl, bl_old, turb, fluid_params, grid)
call add_coriolis(bl, ug, vg, f_cor, dt, grid)
time_current = time_current + dt
if (time_current >=out_dt) then
......@@ -237,11 +303,12 @@ 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 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) = umst
series_val(3) = turb%ustr
series_val(4) = hpbl
series_val(5) = 0.0
......@@ -275,10 +342,9 @@ program gabls1
call close_file(scan_f)
call deallocate_pbl_grid(grid)
!deallocate(ttemp)
!deallocate(ttime)
!deallocate(utemp)
!deallocate(vtemp)
call scm_data_deallocate(bl)
call scm_data_deallocate(bl_old)
call pbl_data_deallocate(turb)
end program gabls1
......@@ -288,7 +354,10 @@ subroutine init_state(bl, ug_,vg_,tref)
implicit none
real(kind=rf):: ug_,vg_,tref
type(stateBLDataType), intent(out) :: bl
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
......@@ -372,7 +441,7 @@ subroutine get_surface_from_config(model, surface, z0,zh )
call c_config_get_string("model.id"//C_NULL_CHAR, config_field, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
stop
end if
model = get_model_id(char_array2str(config_field))
if (model == -1) then
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment