added solver and turb functions

CMakeLists.txt

@@ -35,6 +35,7 @@ set(lib_files
     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})

src/pbl_fo_turb.f90

+! 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

src/pbl_solver.f90

+! 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

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

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