From ed08e8fc3b6bb77a60a46d31a3b87f6d65142f12 Mon Sep 17 00:00:00 2001 From: Andrey Debolskiy <and.debol@gmail.com> Date: Fri, 14 Feb 2025 06:46:08 +0300 Subject: [PATCH] reworked some functions --- CMakeLists.txt | 9 +- config-examples/config-gabls1-ex.txt | 34 +++++-- src/config-utils.f90 | 64 +++++++++++++ src/fo_stability_functions.f90 | 91 +++++++++++++++++++ src/pbl_fo_turb.f90 | 121 +++++++++++++++++++++---- src/pbl_grid.f90 | 1 + src/pbl_turb_common.f90 | 112 +++++++++++++++++++++++ src/pbl_turb_data.f90 | 52 ++++++++--- src/scm_state_data.f90 | 49 ++++++---- src/state_utils.f90 | 67 +++++++++++++- src/tests/gabls1.f90 | 131 ++++++++++++++++++++------- 11 files changed, 640 insertions(+), 91 deletions(-) create mode 100644 src/fo_stability_functions.f90 create mode 100644 src/pbl_turb_common.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index a6c7d11..46f8633 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -35,7 +35,10 @@ set(lib_files src/state_utils.f90 src/scm_state_data.f90 src/pbl_turb_data.f90 - src/pbl_solver.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() \ No newline at end of file + +endif() + diff --git a/config-examples/config-gabls1-ex.txt b/config-examples/config-gabls1-ex.txt index 5f88583..5cf74b0 100644 --- a/config-examples/config-gabls1-ex.txt +++ b/config-examples/config-gabls1-ex.txt @@ -1,12 +1,16 @@ 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; - } +} fluid { pref= 1013.4; #hPa @@ -14,16 +18,28 @@ fluid { g = 9.81; # m/s2 beta = g /tref; kappa = 0.4; # von Karman constant - } +} -surface_model { +grid { + type = "uniform"; - id = "esm"; # optional: "esm" (default), "sheba", "most", "log" - type = "land"; - z0_m = 0.01; # aerodynamic roughness [m] - z0_h = -1; # no prescribed value # -> using scheme assigned by surface type - } + 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] + 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 diff --git a/src/config-utils.f90 b/src/config-utils.f90 index d302626..0076d75 100644 --- a/src/config-utils.f90 +++ b/src/config-utils.f90 @@ -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) diff --git a/src/fo_stability_functions.f90 b/src/fo_stability_functions.f90 new file mode 100644 index 0000000..0f32177 --- /dev/null +++ b/src/fo_stability_functions.f90 @@ -0,0 +1,91 @@ +! 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 diff --git a/src/pbl_fo_turb.f90 b/src/pbl_fo_turb.f90 index 959dbcc..cf617e1 100644 --- a/src/pbl_fo_turb.f90 +++ b/src/pbl_fo_turb.f90 @@ -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 diff --git a/src/pbl_grid.f90 b/src/pbl_grid.f90 index 90a0d0a..cb1fe27 100644 --- a/src/pbl_grid.f90 +++ b/src/pbl_grid.f90 @@ -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)) diff --git a/src/pbl_turb_common.f90 b/src/pbl_turb_common.f90 new file mode 100644 index 0000000..948abbc --- /dev/null +++ b/src/pbl_turb_common.f90 @@ -0,0 +1,112 @@ +! 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 diff --git a/src/pbl_turb_data.f90 b/src/pbl_turb_data.f90 index 9e55a8a..7f6e88d 100644 --- a/src/pbl_turb_data.f90 +++ b/src/pbl_turb_data.f90 @@ -1,16 +1,43 @@ ! 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 diff --git a/src/scm_state_data.f90 b/src/scm_state_data.f90 index ce08dc9..bf65975 100644 --- a/src/scm_state_data.f90 +++ b/src/scm_state_data.f90 @@ -1,21 +1,31 @@ ! 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 diff --git a/src/state_utils.f90 b/src/state_utils.f90 index d4c3af3..8c19a71 100644 --- a/src/state_utils.f90 +++ b/src/state_utils.f90 @@ -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 diff --git a/src/tests/gabls1.f90 b/src/tests/gabls1.f90 index 0ddb027..d4b852a 100644 --- a/src/tests/gabls1.f90 +++ b/src/tests/gabls1.f90 @@ -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,42 +146,88 @@ 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) + 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 - 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 -- GitLab