416 lines
7.8 KiB
Fortran
416 lines
7.8 KiB
Fortran
module m_terms
|
|
|
|
use m_parameters
|
|
use m_arrays
|
|
use m_calculate
|
|
|
|
implicit none
|
|
|
|
real*8, allocatable, dimension(:,:,:) :: ddx_nx
|
|
real*8, allocatable, dimension(:,:,:) :: t1_fsd
|
|
real*8, allocatable, dimension(:,:,:) :: t2_fsd
|
|
real*8, allocatable, dimension(:,:,:) :: ddz_c
|
|
real*8, allocatable, dimension(:,:,:) :: vn
|
|
real*8, allocatable, dimension(:,:,:) :: ddy_c
|
|
real*8, allocatable, dimension(:,:,:) :: ddx_c
|
|
real*8, allocatable, dimension(:) :: fsd_avg_fu
|
|
real*8, allocatable, dimension(:) :: fsd_avg_t2
|
|
real*8, allocatable, dimension(:) :: fsd_avg_t1
|
|
real*8, allocatable, dimension(:,:,:) :: ny
|
|
real*8, allocatable, dimension(:,:,:) :: wrate
|
|
real*8, allocatable, dimension(:) :: avg_fsd
|
|
real*8, allocatable, dimension(:) :: fsd_avg_absk
|
|
real*8, allocatable, dimension(:,:,:) :: nx
|
|
real*8, allocatable, dimension(:,:,:) :: ddy_ny
|
|
real*8, allocatable, dimension(:,:,:) :: nz
|
|
real*8, allocatable, dimension(:,:,:) :: ddz_nz
|
|
real*8, allocatable, dimension(:,:,:) :: fu
|
|
real*8, allocatable, dimension(:,:,:) :: c
|
|
real*8, allocatable, dimension(:,:,:) :: fsd
|
|
real*8, allocatable, dimension(:,:,:) :: d2dy_c
|
|
real*8, allocatable, dimension(:,:,:) :: d2dx_c
|
|
real*8, allocatable, dimension(:,:,:) :: d2dz_c
|
|
real*8, allocatable, dimension(:,:,:) :: absk
|
|
real*8, allocatable, dimension(:,:,:) :: sd
|
|
|
|
|
|
contains
|
|
|
|
subroutine m_terms_init
|
|
|
|
integer :: ierr
|
|
|
|
allocate(ddx_nx(nxp,nyp,nzp), stat=ierr) ; ddx_nx = 0.
|
|
allocate(t1_fsd(nxp,nyp,nzp), stat=ierr) ; t1_fsd = 0.
|
|
allocate(t2_fsd(nxp,nyp,nzp), stat=ierr) ; t2_fsd = 0.
|
|
allocate(ddz_c(nxp,nyp,nzp), stat=ierr) ; ddz_c = 0.
|
|
allocate(vn(nxp,nyp,nzp), stat=ierr) ; vn = 0.
|
|
allocate(ddy_c(nxp,nyp,nzp), stat=ierr) ; ddy_c = 0.
|
|
allocate(ddx_c(nxp,nyp,nzp), stat=ierr) ; ddx_c = 0.
|
|
allocate(fsd_avg_fu(nxp), stat=ierr) ; fsd_avg_fu = 0.
|
|
allocate(fsd_avg_t2(nxp), stat=ierr) ; fsd_avg_t2 = 0.
|
|
allocate(fsd_avg_t1(nxp), stat=ierr) ; fsd_avg_t1 = 0.
|
|
allocate(ny(nxp,nyp,nzp), stat=ierr) ; ny = 0.
|
|
allocate(wrate(nxp,nyp,nzp), stat=ierr) ; wrate = 0.
|
|
allocate(avg_fsd(nxp), stat=ierr) ; avg_fsd = 0.
|
|
allocate(fsd_avg_absk(nxp), stat=ierr) ; fsd_avg_absk = 0.
|
|
allocate(nx(nxp,nyp,nzp), stat=ierr) ; nx = 0.
|
|
allocate(ddy_ny(nxp,nyp,nzp), stat=ierr) ; ddy_ny = 0.
|
|
allocate(nz(nxp,nyp,nzp), stat=ierr) ; nz = 0.
|
|
allocate(ddz_nz(nxp,nyp,nzp), stat=ierr) ; ddz_nz = 0.
|
|
allocate(fu(nxp,nyp,nzp), stat=ierr) ; fu = 0.
|
|
allocate(c(nxp,nyp,nzp), stat=ierr) ; c = 0.
|
|
allocate(fsd(nxp,nyp,nzp), stat=ierr) ; fsd = 0.
|
|
allocate(d2dy_c(nxp,nyp,nzp), stat=ierr) ; d2dy_c = 0.
|
|
allocate(d2dx_c(nxp,nyp,nzp), stat=ierr) ; d2dx_c = 0.
|
|
allocate(d2dz_c(nxp,nyp,nzp), stat=ierr) ; d2dz_c = 0.
|
|
allocate(absk(nxp,nyp,nzp), stat=ierr) ; absk = 0.
|
|
allocate(sd(nxp,nyp,nzp), stat=ierr) ; sd = 0.
|
|
|
|
|
|
end subroutine m_terms_init
|
|
|
|
|
|
subroutine m_terms_finalize
|
|
|
|
deallocate(ddx_nx)
|
|
deallocate(t1_fsd)
|
|
deallocate(t2_fsd)
|
|
deallocate(ddz_c)
|
|
deallocate(vn)
|
|
deallocate(ddy_c)
|
|
deallocate(ddx_c)
|
|
deallocate(fsd_avg_fu)
|
|
deallocate(fsd_avg_t2)
|
|
deallocate(fsd_avg_t1)
|
|
deallocate(ny)
|
|
deallocate(wrate)
|
|
deallocate(avg_fsd)
|
|
deallocate(fsd_avg_absk)
|
|
deallocate(nx)
|
|
deallocate(ddy_ny)
|
|
deallocate(nz)
|
|
deallocate(ddz_nz)
|
|
deallocate(fu)
|
|
deallocate(c)
|
|
deallocate(fsd)
|
|
deallocate(d2dy_c)
|
|
deallocate(d2dx_c)
|
|
deallocate(d2dz_c)
|
|
deallocate(absk)
|
|
deallocate(sd)
|
|
|
|
|
|
end subroutine m_terms_finalize
|
|
|
|
|
|
subroutine m_terms_calculate_pass1
|
|
|
|
integer :: i, j, k
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
c(i,j,k) = ( 1.0 - y(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call ddz ( ddz_c, c )
|
|
call ddy ( ddy_c, c )
|
|
call d2dx ( d2dx_c, c )
|
|
call ddx ( ddx_c, c )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
fsd(i,j,k) = ( sqrt ( ( ( ((ddx_c(i,j,k))*(ddx_c(i,j,k))) + ((ddy_c(i,j,k))*(ddy_c(i,j,k))) ) + ((ddz_c(i,j,k))*(ddz_c(i,j,k))) ) ) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
wrate(i,j,k) = ( rxn_rate ( c(i,j,k) ) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call d2dy ( d2dy_c, c )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
avg_fsd(i) = avg_fsd(i) + fsd(i,j,k)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
nx(i,j,k) = ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call ddx ( ddx_nx, nx )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
ny(i,j,k) = ( ( - ddy_c(i,j,k) ) / fsd(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call ddy ( ddy_ny, ny )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
nz(i,j,k) = ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
vn(i,j,k) = ( ( ( u(i,j,k) * nx(i,j,k) ) + ( v(i,j,k) * ny(i,j,k) ) ) + ( w(i,j,k) * nz(i,j,k) ) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call ddz ( ddz_nz, nz )
|
|
call d2dz ( d2dz_c, c )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
sd(i,j,k) = ( ( ( ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) ) * rod ) + wrate(i,j,k) ) / fsd(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
fu(i,j,k) = ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
fsd_avg_fu(i) = fsd_avg_fu(i) + fu(i,j,k) * fsd(i,j,k)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
absk(i,j,k) = ( dabs ( ( ( ddx_nx(i,j,k) + ddy_ny(i,j,k) ) + ddz_nz(i,j,k) ) ) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
fsd_avg_absk(i) = fsd_avg_absk(i) + absk(i,j,k) * fsd(i,j,k)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
|
|
end subroutine m_terms_calculate_pass1
|
|
|
|
|
|
subroutine m_terms_average_pass1 (nfiles)
|
|
|
|
integer :: nfiles
|
|
real*8 :: denum
|
|
|
|
denum=real(nfiles*nyp*nzp)
|
|
|
|
avg_fsd = avg_fsd / denum
|
|
fsd_avg_fu = fsd_avg_fu / denum / avg_fsd
|
|
fsd_avg_absk = fsd_avg_absk / denum / avg_fsd
|
|
|
|
|
|
end subroutine m_terms_average_pass1
|
|
|
|
|
|
subroutine m_terms_calculate_pass2
|
|
|
|
integer :: i, j, k
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
c(i,j,k) = ( 1.0 - y(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call ddz ( ddz_c, c )
|
|
call ddy ( ddy_c, c )
|
|
call ddx ( ddx_c, c )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
fsd(i,j,k) = ( sqrt ( ( ( ((ddx_c(i,j,k))*(ddx_c(i,j,k))) + ((ddy_c(i,j,k))*(ddy_c(i,j,k))) ) + ((ddz_c(i,j,k))*(ddz_c(i,j,k))) ) ) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
ny(i,j,k) = ( ( - ddy_c(i,j,k) ) / fsd(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
wrate(i,j,k) = ( rxn_rate ( c(i,j,k) ) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call d2dy ( d2dy_c, c )
|
|
call d2dx ( d2dx_c, c )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
nx(i,j,k) = ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call ddx ( ddx_nx, nx )
|
|
call ddy ( ddy_ny, ny )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
nz(i,j,k) = ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
vn(i,j,k) = ( ( ( u(i,j,k) * nx(i,j,k) ) + ( v(i,j,k) * ny(i,j,k) ) ) + ( w(i,j,k) * nz(i,j,k) ) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
call ddz ( ddz_nz, nz )
|
|
call d2dz ( d2dz_c, c )
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
absk(i,j,k) = ( dabs ( ( ( ddx_nx(i,j,k) + ddy_ny(i,j,k) ) + ddz_nz(i,j,k) ) ) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
t2_fsd(i,j,k) = (((absk(i,j,k) - fsd_avg_absk(i)))*((absk(i,j,k) - fsd_avg_absk(i))))
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
fsd_avg_t2(i) = fsd_avg_t2(i) + t2_fsd(i,j,k) * fsd(i,j,k)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
sd(i,j,k) = ( ( ( ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) ) * rod ) + wrate(i,j,k) ) / fsd(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
fu(i,j,k) = ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
t1_fsd(i,j,k) = (((fu(i,j,k) - fsd_avg_fu(i)))*((fu(i,j,k) - fsd_avg_fu(i))))
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do k = 1, nzp
|
|
do j = 1, nyp
|
|
do i = 1, nxp
|
|
fsd_avg_t1(i) = fsd_avg_t1(i) + t1_fsd(i,j,k) * fsd(i,j,k)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
|
|
end subroutine m_terms_calculate_pass2
|
|
|
|
|
|
subroutine m_terms_average_pass2 (nfiles)
|
|
|
|
integer :: nfiles
|
|
real*8 :: denum
|
|
|
|
denum=real(nfiles*nyp*nzp)
|
|
|
|
fsd_avg_t1 = fsd_avg_t1 / denum / avg_fsd
|
|
fsd_avg_t2 = fsd_avg_t2 / denum / avg_fsd
|
|
|
|
|
|
end subroutine m_terms_average_pass2
|
|
|
|
|
|
end module m_terms
|
|
|