Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
! 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
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:kbltop-1) = 0.0
bb(1:kbltop-1) = 0.0
cc(1:kbltop-1) = 0.0
k = kbltop
dtdz = dt / (grid%dzc(k))
aa(k) = 0
cc(k) = (kdiff(k)/rho(k)) * dtdz / grid%dze(k)
write(*,*) 'here_diff1.5', kdiff(k), kbltop
write(*,*) 'KTVDM', k, aa(k), bb(k), cc(k), rho(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)
!bottom boundary
k = grid%kmax
dtdz = dt / (grid%dzc(k))
aa(k) = (kdiff(k-1)/rho(k)) * dtdz / grid%dze(k-1)
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
real, allocatable:: aa(:), bb(:), cc(:), ff(:)
real, allocatable:: prgna(:), prgnz(:)
integer k, integer, ktop, kmax
kmax = grid%kmax
write(*,*) 'here_diff'
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
write(*,*) 'here_diff1', bl%kpbl
ktop = bl%kpbl
do k = bl%kpbl-1,1,-1
if(bl%vdcuv(k) > 0.e0) then
ktop = k
end if
enddo
! fill for temperature and specific humidity
call fill_tridiag(aa, bb, cc, bl%rho, bl%vdctq, ktop, 0.0, grid, dt)
write(*,*) 'here_diff2'
call factorize_tridiag(ktop, kmax, aa, bb, cc, prgna, prgnz)
do k = ktop,kmax-1
ff(k) = bl%theta(k)
write(*,*) '2', k, aa(k), bb(k), cc(k), ff(k), grid%dze(k)
end do
ff(kmax) = bl%theta(kmax) + (dt/kmax) * bl%surf%hs /bl%rho(kmax)
write(*,*) 'here_diff3'
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%theta)
write(*,*) 'here_diff4'
do k = ktop,kmax-1
ff(k) = bl%qv(k)
end do
ff(kmax) = bl%qv(kmax) + (dt/kmax) * bl%surf%es /bl%rho(kmax)
call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%qv)
write(*,*) 'here_diff4'
!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)
end subroutine solve_diffusion
end module pbl_solver