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