diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5df81955e23cc0dd653e531e852d5aee9a188ce1..a6c7d111531d8ed72ac61a73a9da762846a22b88 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -31,10 +31,11 @@ endif(BUILD_DOC)
set(lib_files
src/parkinds.f90
src/pbl_grid.f90
- src/phys_fluid.f90
+ src/phys_fluid.f90
src/state_utils.f90
src/scm_state_data.f90
src/pbl_turb_data.f90
+ src/pbl_solver.f90
src/diag_pbl.f90)
add_library(abl_lib ${lib_files})
diff --git a/src/pbl_fo_turb.f90 b/src/pbl_fo_turb.f90
new file mode 100644
index 0000000000000000000000000000000000000000..fe140e683046c6692f05c15e51f2641a0f4e8e3d
--- /dev/null
+++ b/src/pbl_fo_turb.f90
@@ -0,0 +1,11 @@
+! Created by Andrey Debolskiy on 29.11.2024.
+
+module pbl_fo_turb
+ implicit none
+
+ public get_fo_coeffs
+ public get_turb_length
+ public
+ contains
+
+end module pbl_fo_turb
\ No newline at end of file
diff --git a/src/pbl_solver.f90 b/src/pbl_solver.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ce9530e9d99e166bf1a043ef63ba98b9de508187
--- /dev/null
+++ b/src/pbl_solver.f90
@@ -0,0 +1,109 @@
+! Created by Andrey Debolskiy on 29.11.2024.
+
+module pbl_solver
+ use parkinds, only: rf=>kind_rf, im=>kind_im
+ use scm_state_data
+ use pbl_turb_data
+ implicit none
+
+ public factorize_tridiag
+ public solve_tridiag
+ public fill_tridiag
+ public solve_diffusion
+ contains
+ subroutine factorize_tridiag(ktvd, kl, aa, bb, cc, prgna, prgnz)
+ implicit none
+ integer, intent(in) :: ktvd, kl
+ real, intent(in), dimension(kl) :: aa, bb, cc
+ real, intent(out), dimension(kl) :: prgna, prgnz
+
+ integer :: k
+
+ prgna(ktvd) = cc(ktvd) / bb(ktvd)
+ do k = ktvd+1, kl
+ prgnz(k) = 1.0e0 / (bb(k) - aa(k) * prgna(k-1))
+ prgna(k) = cc(k) * prgnz(k)
+ end do
+ end subroutine factorize_tridiag
+
+ !> reduce tridiagonal system to bidiagonal after matrix factorization
+ !! - bb(k) y(k) + cc(k) y(k+1) = -f(k), k = ktvd
+ !! aa(k) y(k-1) - bb(k) y(k) + cc(k) y(k+1) = -f(k), k = ktvd+1..kl-1
+ !! aa(k) y(k-1) - bb(k) y(k) = -f(k), k = kl
+ !! assuming cc(kl) = 0.0
+ !! reduced system is
+ !! y(k) - prgna(k) y(k+1) = prgnb(k), k = ktvd..kl-1
+ !! y(k) = prgnb(k), k = kl
+ !! then solve via backward substitution
+ subroutine solve_tridiag(ktvd, kl, aa, bb, cc, ff, prgna, prgnz, y)
+ implicit none
+ integer, intent(in) :: ktvd, kl
+ real, intent(in), dimension(kl) :: aa, bb, cc, ff
+ real, intent(in), dimension(kl) :: prgna, prgnz
+ real, intent(out), dimension(kl) :: y
+
+ integer :: k
+ real :: prgnb(kl)
+
+ prgnb(ktvd) = ff(ktvd) / bb(ktvd)
+ do k = ktvd+1, kl
+ prgnb(k) = prgnz(k) * (ff(k) + aa(k) * prgnb(k-1))
+ end do
+ y(kl) = prgnb(kl)
+ do k = kl-1, ktvd, -1
+ y(k) = prgna(k) * y(k+1) + prgnb(k)
+ end do
+ end subroutine solve_tridiag
+
+ subroutine fill_tridiag(aa, bb, cc, rho, kdiff, kbltop, grid, dt)
+ use pbl_grid, only: pblgridDataType
+ implicit none
+ real, dimension(*), intent(in):: rho, kdiff
+ real, intent(in):: dt
+ integer, intent(in):: kbltop
+ type(pblgridDataType),intent(in):: grid
+
+ real, dimension(*), intent(out):: aa, bb, cc
+
+ real:: dtdz
+ integer:: k
+
+ !nulify before top boundary
+ aa(1:kbltop) = 0.0
+ bb(1:kbltop) = 0.0
+ cc(1:kbltop) = 0.0
+ !top boundary condition: flux = 0
+ k = kbltop
+ dtdz = dt / (grid%dzc(k))
+ aa(k) = 0
+ cc(k) = (kdiff(k)/rho(k)) * dtdz / grid%dze(k)
+
+ do k = kbltop + 1, grid%kmax -1
+ dtdz = dt / (grid%dzc(k))
+ aa(k) = (kdiff(k - 1)/rho(k)) * dtdz / grid%dze(k-1)
+ cc(k) = (kdiff(k)/rho(k)) * dtdz / grid%dze(k)
+ bb(k) = 1.0 + aa(k) + cc(k)
+ end do
+
+ !bottom boundary
+ k = grid%kmax
+ dtdz = dt / (grid%dzc(k))
+ aa(k) = (kdiff(k-1)/rho(k)) * dtdz / grid%dze(k-1)
+ bb(k) = 1.0 + aa(k)
+ cc(k) = 0.0
+ end subroutine fill_tridiag
+
+ subroutine solve_diffusion(bl, bl_old, turb, fluid, grid)
+ use scm_state_data, only : stateBLDataType
+ use pbl_turb_data, only : turbBLDataType
+ use phys_fluid, only: fluidParamsDataType
+ use pbl_grid, only : pblgridDataType
+ implicit none
+ type(stateBLDataType), intent(out):: bl
+ type(stateBLDataType), intent(in):: bl_old
+ type(turbBLDataType), intent(in):: turb
+ type(fluidParamsDataType), intent(in) :: fluid
+ type(pblgridDataType), intent(in) :: grid
+
+ end subroutine solve_diffusion
+end module pbl_solver
\ No newline at end of file
diff --git a/src/state_utils.f90 b/src/state_utils.f90
index 311241cc30eb23d4c4cba32bcc00c3605cba2f2b..2d0eebcfe5e874660856f6a1d28e3b6fc7abf9c4 100644
--- a/src/state_utils.f90
+++ b/src/state_utils.f90
@@ -2,7 +2,11 @@
module state_utils
-
+ public geo
+ public theta2ta
+ public ta2theta
+ public dqsdt_sc, qmax_sc, qmax
+ public DQSDT2, DQSDT, CLDT
contains
SUBROUTINE GEO(T,FS,F,SIG, R, LL)
@@ -31,23 +35,6 @@ module state_utils
-
- SUBROUTINE NANOLOVKA(X,INAN)
- implicit none
- real, intent(in):: x
- integer, intent(out):: inan
- inan=1
-
- IF(X.LE.0.0) THEN
- INAN=0
- END IF
-
- IF(X.GE.0.0) THEN
- INAN=0
- END IF
-
- END SUBROUTINE NANOLOVKA
-
subroutine theta2ta(ta, theta, p0, sig, appa, kl)
implicit none
real, dimension(kl), intent(in):: theta
@@ -77,4 +64,76 @@ module state_utils
theta(k) = ta(k) * (p0 * sig(k))**(-appa)
end do
end subroutine ta2theta
+
+ real function qmax_sc(t, p, aa)
+ implicit none
+
+ real, intent(in) :: t, p, aa
+ real :: pb2, pb
+
+ pb2 = t - 273.2e0
+ if ((aa.gt.1.0e0.and.pb2.ge.0.0e0).or.(aa.lt.1.0e0.and.pb2.ge.-12.0e0)) then
+ pb=17.55e0/(t-31.1e0)
+ else
+ pb=21.85e0/(t-7.65e0)
+ end if
+ qmax_sc=3.80042e-3*exp(pb*pb2)/p
+ end function qmax_sc
+ real function qmax(T,P, AA )
+ IMPLICIT NONE
+ real, intent(in):: t, p, aa
+ !local
+ real pb2, pb
+ PB2=T-273.2E0
+ IF((AA.GT.1.0E0.AND.PB2.GE.0.0E0).OR. &
+ (AA.LT.1.0E0.AND.PB2.GE.-12.0E0))THEN
+ PB=17.55E0/(T-31.1E0 )
+ ELSE
+ PB=21.85E0/(T-7.65E0)
+ ENDIF
+ QMAX=3.80042E-3*EXP(PB*PB2)/P
+ RETURN
+ end function qmax
+
+
+ real function dqsdt_sc(a, b, h)
+ real, intent(in) :: a, b, h
+
+ dqsdt_sc=4248.855e0*h*qmax_sc(a,b, 0.0e0)/(a-31.10e0)**2
+ end function dqsdt_sc
+
+ REAL FUNCTION CLDT(Q,T,P,C)
+ !
+ IMPLICIT NONE
+ REAL, INTENT(IN) :: Q,T,P,C
+ ! REAL DQSDT
+ !REAL :: QMAX
+
+ CLDT=(Q-C*QMAX(T,P, 0.0E0))/(1.0E0+2488.851E0*DQSDT(T,P,C))
+ RETURN
+ END FUNCTION CLDT
+
+ REAL FUNCTION DQSDT(A,B,H)
+ !
+ IMPLICIT NONE
+ REAL, INTENT(IN) :: A, B, H
+ !REAL :: QMAX
+
+ DQSDT=4248.855E0*H*QMAX(A,B, 0.0E0)/(A-31.10E0)**2
+ RETURN
+ END FUNCTION DQSDT
+
+ REAL FUNCTION DQSDT2(A,B,H)
+
+ IMPLICIT NONE
+
+ REAL, INTENT(IN) :: A, B, H
+ ! REAL DQSDT
+ !REAL :: QMAX
+
+ DQSDT2=4248.855E0*H*(DQSDT(A,B, H )/(A-31.10E0)**2 - &
+ 2.0E0*QMAX (A,B, 0.0E0)/(A-31.10E0)**3)
+ RETURN
+ END FUNCTION DQSDT2
+
end module state_utils
\ No newline at end of file
diff --git a/src/tests/gabls1.f90 b/src/tests/gabls1.f90
index 91478bce2a6c83b0676d29679f4967e6d81c7c8a..8bd572b69dec5a21ddfd4604cd6ffa9ae1cc4319 100644
--- a/src/tests/gabls1.f90
+++ b/src/tests/gabls1.f90
@@ -10,6 +10,7 @@ program gabls1
use scm_sfx_data, only: taux, tauy, umst, hflx, qflx, cu
use scm_state_data
use pbl_turb_data
+ use pbl_solver
use state_utils, only : geo, theta2ta, ta2theta
use diag_pbl
use pbl_grid
@@ -38,6 +39,7 @@ program gabls1
type(pblgridDataType) grid
type(stateBLDataType):: bl, bl_old;
+ type(turbBLDataType):: turb
type(meteoDataType) :: meteo_cell
type(sfxDataType) :: sfx_cell
@@ -191,6 +193,8 @@ program gabls1
END DO
call theta2ta(bl_old%temp, bl_old%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
+
+
!io setup
nt = floor((time_end - time_begin)/output_dt)
allocate(ttime(nt))
@@ -217,7 +221,8 @@ program gabls1
hflx = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U * meteo_cell%dT
umst = sfx_cell%Cm * meteo_cell%U
- !call vdiff_scm()
+ 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