Newer
Older
! 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
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, cm2u, grid, dt)
use pbl_grid, only: pblgridDataType
implicit none
real, dimension(*), intent(in):: rho, kdiff
integer, intent(in):: kbltop
type(pblgridDataType),intent(in):: grid
real:: dtdz
integer:: k
!nulify before top boundary
aa(1:grid%kmax) = 0.0
bb(1:grid%kmax) = 0.0
cc(1:grid%kmax) = 0.0
!top boundary condition: flux = 0
k = kbltop
dtdz = dt / (grid%z_edge(k) - grid%z_edge(k-1) )
cc(k) = (kdiff(k)/rho(k)) * dtdz / (grid%z_cell(k+1) - grid%z_cell(k))
dtdz = dt / (grid%z_edge(k) - grid%z_edge(k-1) )
aa(k) = (kdiff(k - 1)/rho(k)) * dtdz / (grid%z_cell(k) - grid%z_cell(k-1))
cc(k) = (kdiff(k)/rho(k)) * dtdz /(grid%z_cell(k+1) - grid%z_cell(k))
end do
!bottom boundary
k = grid%kmax
dtdz = dt / (grid%z_edge(k) - grid%z_edge(k-1) )
aa(k) = (kdiff(k-1)/rho(k)) * dtdz / (grid%z_cell(k) - grid%z_cell(k-1))
bb(k) = 1.0 + aa(k) - dtdz * cm2u / rho(k)
cc(k) = 0.0
end subroutine fill_tridiag
subroutine solve_diffusion(bl, bl_old, turb, 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(stateBLDataType), intent(in):: bl_old
type(turbBLDataType), intent(in):: turb
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
real, save, allocatable:: aa(:), bb(:), cc(:), ff(:)
real, allocatable:: prgna(:), prgnz(:)
integer k, integer, ktop, kmax
kmax = grid%kmax
if (.not.(allocated(aa))) then
allocate(aa(grid%kmax), source=0.0)
end if
if (.not.(allocated(bb))) then
allocate(bb(grid%kmax), source=0.0)
end if
if (.not.(allocated(cc))) then
allocate(cc(grid%kmax), source=0.0)
end if
if (.not.(allocated(ff))) then
allocate(ff(grid%kmax), source=0.0)
end if
if (.not.(allocated(prgna))) then
allocate(prgna(grid%kmax), source=0.0)
end if
if (.not.(allocated(prgnz))) then
allocate(prgnz(grid%kmax), source=0.0)
end if
ktop = bl%kpbl
!this loop is generally not needed
!do k = bl%kpbl-1,1,-1
! if(bl%vdcuv(k) > 0.e0) then
! ktop = k
! end if
!enddo
call fill_tridiag(aa, bb, cc, bl%rho, bl%vdctq, ktop, 0.0, grid, dt)
call factorize_tridiag(ktop, kmax, aa, bb, cc, prgna, prgnz)
do k = ktop,kmax-1
ff(k) = bl%theta(k)
end do
ff(kmax) = bl%theta(kmax) &
+ (dt/grid%dze(kmax)) * bl%surf%hs / (bl%rho(kmax))
write(*,*) 'ff'!, ff!, aa,bb,cc
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%theta)
ff(kmax) = bl%qv(kmax) + (dt/grid%dze(kmax)) * bl%surf%es /bl%rho(kmax)
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%qv)
!velocity
call fill_tridiag(aa, bb, cc, bl%rho, bl%vdcuv, ktop, bl%surf%cm2u, grid, dt)
call factorize_tridiag(ktop, kmax, aa, bb, cc, prgna, prgnz)
do k = ktop,kmax-1
ff(k) = bl%u(k)
end do
ff(kmax) = bl%u(kmax)
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%u)
do k = ktop,kmax-1
ff(k) = bl%v(k)
end do
ff(kmax) = bl%v(kmax)
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%v)
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
subroutine apply_subsidence(bl, w_s, 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(stateBLDataType), intent(inout):: bl
type(fluidParamsDataType), intent(in) :: fluid
type(pblgridDataType), intent(in) :: grid
real, intent(in):: dt
real, dimension(*), intent(in):: w_s
real, save, allocatable :: subs_u(:), subs_v(:), &
subs_theta(:), subs_qv(:)
real:: dtdz, rhom, rhop
integer k, kmax
kmax = grid%kmax
if (.not.(allocated(subs_u))) then
allocate(subs_u(grid%kmax), source=0.0)
end if
if (.not.(allocated(subs_v))) then
allocate(subs_v(grid%kmax), source=0.0)
end if
if (.not.(allocated(subs_theta))) then
allocate(subs_theta(grid%kmax), source=0.0)
end if
if (.not.(allocated(subs_qv))) then
allocate(subs_qv(grid%kmax), source=0.0)
end if
! central differences
do k = kmax-1, 2, -1
rhop = 0.5*(bl%rho(k) + bl%rho(k+1))
rhom = 0.5*(bl%rho(k) + bl%rho(k-1))
subs_u(k) = (0.5/bl%rho(k)) &
* (rhop * w_s(k) * (bl%u(k)-bl%u(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1)) &
+ &
rhom * w_s(k-1) * (bl%u(k)-bl%u(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1)))
subs_v(k) = (0.5/bl%rho(k)) &
* (rhop * w_s(k) * (bl%v(k)-bl%v(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1)) &
+ &
rhom * w_s(k-1) * (bl%v(k)-bl%v(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1)))
subs_theta(k) = (0.5/bl%rho(k)) &
* (rhop * w_s(k) * (bl%theta(k)-bl%theta(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1)) &
+ &
rhom * w_s(k-1) * (bl%theta(k)-bl%theta(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1)))
subs_qv(k) = (0.5/bl%rho(k)) &
* (rhop * w_s(k) * (bl%qv(k)-bl%qv(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1)) &
+ &
rhom * w_s(k-1) * (bl%qv(k)-bl%qv(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1)))
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))
subs_v(k) = w_s(k-1) * (bl%v(k)-bl%v(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1))
subs_theta(k) = w_s(k-1) * (bl%theta(k)-bl%theta(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1))
subs_qv(k) = w_s(k-1) * (bl%qv(k)-bl%qv(k-1)) &
/(grid%z_cell(k)-grid%z_cell(k-1))
! top boundary
k = 1
subs_u(k) = w_s(k+1) * (bl%u(k)-bl%u(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1))
subs_v(k) = w_s(k+1) * (bl%v(k)-bl%v(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1))
subs_theta(k) = w_s(k+1) * (bl%theta(k)-bl%theta(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1))
subs_qv(k) = w_s(k+1) * (bl%qv(k)-bl%qv(k+1)) &
/(grid%z_cell(k)-grid%z_cell(k+1))
!apply
do k=kmax, 1, -1
bl%u(k) = bl%u(k) - subs_u(k) * dt
bl%v(k) = bl%v(k) - subs_v(k) * dt
bl%theta(k) = bl%theta(k) - subs_theta(k) * dt
bl%qv(k) = bl%qv(k) - subs_qv(k) * dt
end do
end subroutine apply_subsidence