把这个计算标准差优化一下,先用集合对背景场作差,然后挑选出1,99%分位数的范围设为有效范围,然后对有效点计算标准差;然后,把那些无效点全部归为0;再往后进行标准化
subroutine normalize_local_ensemble(local_back, local_ensemble, &
i_grid, j_grid, layer_std)
use module_grid
implicit none
real, dimension(:,:,:) :: local_back
real, dimension(:,:,:,:) :: local_ensemble
integer, intent(in) :: i_grid, j_grid
real, dimension(:) :: layer_std
real :: std_temp
integer :: i, j, k, member, istart, iend, jstart, jend, kstart, kend
real :: local_number_of_points, vensemble, vback, delta_en_ba
istart = 1
iend = 2*local_radius+1
jstart = 1
jend = 2*local_radius+1
kstart = 1
kend = kme-kms+1
! 为每个垂直层计算标准差
! do k = kstart, kend
! layer_std(k) = 0.0
! end do
layer_std = 0.0
! print*, 'x from', istart, 'to', iend
! print*, 'y from', jstart, 'to', jend
! 计算每个垂直层的标准差
! write(*, '(A,I3,A,I3)') 'Calculating standard deviation for grid point:', i_grid, ' ', j_grid
! allocate(layer_count(kms:kme))
! do k = kstart, kend
! layer_count(k) = 0.0
! end do
! 向量化计算标准差 - 严格按循环顺序:i-k-j-member
do i = istart, iend
do k = kstart, kend
do j = jstart, jend
do member = 1, nsample
vensemble = local_ensemble(i, k, j, member)
vback = local_back(i, k, j)
! print*, 'ijk=', i, j, k, ' member=', member, &
! ' vensemble=', vensemble, ' vback=', vback
delta_en_ba = vensemble - vback
! if (abs(delta_en_ba) > 30.0) then
! print*, vback, ' vback'
! print*, vensemble, ' vensemble'
! end if
! if (vensemble >= vmin .and. vensemble <= vmax .and. vback >= vmin .and. vback <= vmax) then
layer_std(k) = layer_std(k) + delta_en_ba**2
! layer_count(k) = layer_count(k) + 1.0
! end if
end do
end do
end do
end do
local_number_of_points = real((2*local_radius+1) * (2*local_radius+1) * nsample)
do k = kstart, kend
! 计算标准差
layer_std(k) = sqrt(layer_std(k) / local_number_of_points)
! if (layer_count(k) > 0) then
! layer_std(k) = sqrt(layer_std(k) / layer_count(k))
! else
! layer_std(k) = 0.0
! end if
! 避免除零
! if (layer_std(k) < 1.0e-10) then
! layer_std(k) = 1.0
! end if
! print*, 'k=', k, ' layer_std(k)=', layer_std(k)
! write(*, '(A,I3,A,F10.4,A)') 'Ratio of valid points for layer', k, ':', 100*layer_count(k) / local_number_of_points,'%'
end do
! deallocate(layer_count)
! 对集合预报进行标准化
! write(*, '(A,I3,A,I3)') 'Normalizing ensemble members at grid point:', i_grid, ' ', j_grid
do member = 1, nsample
local_ensemble(istart:iend, kstart:kend, jstart:jend, member) = &
local_ensemble(istart:iend, kstart:kend, jstart:jend, member) - &
local_back(istart:iend, kstart:kend, jstart:jend)
end do
do k = kstart, kend
do member = 1, nsample
if (layer_std(k) == 0.0) then
std_temp = 1.0 ! 避免除零
else
std_temp = layer_std(k)
end if
local_ensemble(istart:iend, k, jstart:jend, member) = &
local_ensemble(istart:iend, k, jstart:jend, member) / std_temp
end do
end do
end subroutine normalize_local_ensemble