Feel free to join the next Helmholtz Hacky Hour #30 on Wednesday, June 16, 2021 from 2PM to 3PM!

Commit d8725316 authored by Stephan Thober's avatar Stephan Thober
Browse files

minor bug changes

- corrected reference date in netCDF output file
- removed print statement for particular cell during debugging
- inverted spatial field to netCDF
parent ee30d840
......@@ -112,15 +112,18 @@ program ED_Kriging
! open netcdf if necessary
call open_netcdf(nc_out, nc_data, nc_time)
do iCell = 1, nCell
! initialize cell
allocate(cell(iCell)%z(jStart:jEnd))
cell(iCell)%z = noDataValue
end do
!$OMP parallel default(shared) &
!$OMP private(iCell)
!$OMP do SCHEDULE(STATIC)
ncellsloop: do iCell = 1, nCell
! initialize cell
allocate(cell(iCell)%z(jStart:jEnd))
cell(iCell)%z = noDataValue
! check DEM
if (nint(cell(iCell)%h) == grid%nodata_value ) then
cell(iCell)%z = gridMeteo%nodata_value
......@@ -145,7 +148,7 @@ program ED_Kriging
do i = 1, gridMeteo%ncols
do j = 1, gridMeteo%nrows
k = k + 1
tmp_array(i, j, :) = cell(k)%z
tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
end do
end do
do i = 1, jEnd - jStart + 1
......
!*************************************************************************
!EmpVar.f90
!AUTHOR: Luis Samaniego-eguiguren
!
!DESCRIPTION
! This program calculates the empirical variogram
! gamma(:,1) : distance
! gamme(:,2) : semivarance
! botm are normalized with gmax(:)
!UPDATES:
! Created: 22.11.2001
! Last update 19.02.2004
!
!**************************************************************************
module mo_EmpVar
implicit none
private
public :: EmpVar
contains
subroutine EmpVar(jd, flagMax)
use mainVar
use mo_kind , only : i4, dp
use kriging
use VarFit
implicit none
integer(i4), intent(in) :: jd
logical, intent(in) :: flagMax
integer(i4) :: i, j, k, nhk, ni
real(dp) :: gk, delta
! ******************************************
! Calculate all combinations (var. function)
! ******************************************
do i=1,nSta-1
do j=i+1,nSta
! clasical
if ( (MetSta(i)%z(jd) /= noDataValue ) .and. &
(MetSta(j)%z(jd) /= noDataValue ) ) then
dz2S(i)%S(j) = (MetSta(i)%z(jd)-MetSta(j)%z(jd))**2
else
dz2S(i)%S(j) = noDataValue
end if
!
! Cressi-Hawkins
!dz2S(i)%S(j)= dsqrt( dabs( MetSta(i)%z(jd) - MetSta(j)%z(jd) ) )
end do
end do
!
! variance h=0 (one pass algorithm)
do i = 1,nSta
if ( MetSta(i)%z(jd) == noDataValue ) cycle
nobs = nobs + 1
delta = MetSta(i)%z(jd) - m0
m0 = m0 + delta / real(nobs, dp)
v0 = v0 + delta * (MetSta(i)%z(jd) - m0)
!write (999,'(4i6, 3f15.4)') y,d,i, nobs, MetSta(i)%z(jd), m0, v0
end do
!
!
! calculate the empiric variogram
!
! selecting pair - bins
if (jd == jStart) then
nbins = ceiling(hMax/dh)
if (.not. allocated (nH)) allocate (Nh(nbins), gamma(nbins,2))
Nh = 0
gamma = 0.0_dp
end if
!
print *, '***WARNING: Removal of outliers in the estimation of the variogram is deactivated'
do i=1,nSta-1
do j=i+1,nSta
if (dz2S(i)%S(j) /= noDataValue ) then
! take values up to max distance
if (dS(i)%S(j) > hMax ) cycle
! ! remove outliers for the estimation of the variogram
! if (flagVarTyp == 2 ) then
! if ( dabs( MetSta(i)%h - MetSta(j)%h ) / dS(i)%S(j) > gradE ) then
! ! write(999,*), 'pair removed', MetSta(i)%id, MetSta(j)%id
! cycle
! end if
! end if
!
k=max(1, ceiling(dS(i)%S(j)/dh))
Nh(k)=Nh(k)+1
gamma(k,2)=gamma(k,2) + dz2S(i)%S(j)
end if
end do
end do
!
! only at the end
if (jd == jEnd) then
!
! consolidate bins N(h)>=30
!forward
do k=1,nbins-1
nhk=Nh(k)
gk=gamma(k,2)
if ((nhk>0) .and. (nhk<30)) then
Nh(k)=-9
gamma(k,2)=0.0_dp
if (Nh(k+1)>0) then
Nh(k+1)=Nh(k+1)+nhk
gamma(k+1,2)=gamma(k+1,2)+gk
else
Nh(k+1)=nhk
gamma(k+1,2)=gk
end if
end if
end do
!backward
do k=nbins,2,-1
nhk=Nh(k)
gk=gamma(k,2)
if ((nhk>0) .and. (nhk<30)) then
Nh(k)=-9
gamma(k,2)=0_dp
if (Nh(k-1)>0) then
Nh(k-1)=Nh(k-1)+nhk
gamma(k-1,2)=gamma(k-1,2)+gk
else
Nh(k-1)=nhk
gamma(k-1,2)=gk
end if
end if
end do
!
ni=0
! estimate semi-variance gamma(i,1)=h, gamma(i,2)=g(h)
do k=1,nbins
if (Nh(k) > 0) then
ni=ni+1
! Classsical
gamma(k,2)=gamma(k,2)/2.0_dp/real(Nh(k), dp)
!
! Cressi-Hawkins: adjust bias
!gamma(k,2)=0.5_dp/(0.457_dp+0.494_dp/dfloat(Nh(k)))*(gamma(k,2)/dfloat(Nh(k)))**4
!
gamma(k,1)=real(k, dp)*dh-real(ni, dp)*dh/2.0_dp
ni=0
else
ni=ni+1
gamma(k,1)=-9.0_dp
end if
end do
!
! scaling
if (flagMax) gmax=maxval(gamma, dim=1)
do k=1,nbins
if (gamma(k,1) > 0.0_dp) then
gamma(k,1) = gamma(k,1)/gmax(1)
gamma(k,2) = gamma(k,2)/gmax(2)
end if
end do
!
!keep variogram
end if
end subroutine EmpVar
end module mo_EmpVar
!*************************************************************************
!EmpVar.f90
!AUTHOR: Luis Samaniego-eguiguren
!
!DESCRIPTION
! This program calculates the empirical variogram
! gamma(:,1) : distance
! gamme(:,2) : semivarance
! botm are normalized with gmax(:)
!UPDATES:
! Created: 22.11.2001
! Last update 19.02.2004
!
!**************************************************************************
module mo_EmpVar
implicit none
private
public :: EmpVar
contains
subroutine EmpVar(jd, flagMax)
use mainVar
use mo_kind , only : i4, dp
use kriging
use VarFit
implicit none
integer(i4), intent(in) :: jd
logical, intent(in) :: flagMax
integer(i4) :: i, j, k, nhk, ni
real(dp) :: gk, delta
! ******************************************
! Calculate all combinations (var. function)
! ******************************************
do i=1,nSta-1
do j=i+1,nSta
! clasical
if ( (MetSta(i)%z(jd) /= noDataValue ) .and. &
(MetSta(j)%z(jd) /= noDataValue ) ) then
dz2S(i)%S(j) = (MetSta(i)%z(jd)-MetSta(j)%z(jd))**2
else
dz2S(i)%S(j) = noDataValue
end if
!
! Cressi-Hawkins
!dz2S(i)%S(j)= dsqrt( dabs( MetSta(i)%z(jd) - MetSta(j)%z(jd) ) )
end do
end do
!
! variance h=0 (one pass algorithm)
do i = 1,nSta
if ( MetSta(i)%z(jd) == noDataValue ) cycle
nobs = nobs + 1
delta = MetSta(i)%z(jd) - m0
m0 = m0 + delta / real(nobs, dp)
v0 = v0 + delta * (MetSta(i)%z(jd) - m0)
!write (999,'(4i6, 3f15.4)') y,d,i, nobs, MetSta(i)%z(jd), m0, v0
end do
!
!
! calculate the empiric variogram
!
! selecting pair - bins
if (jd == jStart) then
nbins = ceiling(hMax/dh)
if (.not. allocated (nH)) allocate (Nh(nbins), gamma(nbins,2))
Nh = 0
gamma = 0.0_dp
end if
!
print *, '***WARNING: Removal of outliers in the estimation of the variogram is activated'
do i=1,nSta-1
do j=i+1,nSta
if (dz2S(i)%S(j) /= noDataValue ) then
! take values up to max distance
if (dS(i)%S(j) > hMax ) cycle
! ! remove outliers for the estimation of the variogram
! if (flagVarTyp == 2 ) then
if ( dabs( MetSta(i)%h - MetSta(j)%h ) / dS(i)%S(j) > gradE ) then
! write(999,*), 'pair removed', MetSta(i)%id, MetSta(j)%id
cycle
end if
! end if
!
k=max(1, ceiling(dS(i)%S(j)/dh))
Nh(k)=Nh(k)+1
gamma(k,2)=gamma(k,2) + dz2S(i)%S(j)
end if
end do
end do
!
! only at the end
if (jd == jEnd) then
!
! consolidate bins N(h)>=30
!forward
do k=1,nbins-1
nhk=Nh(k)
gk=gamma(k,2)
if ((nhk>0) .and. (nhk<30)) then
Nh(k)=-9
gamma(k,2)=0.0_dp
if (Nh(k+1)>0) then
Nh(k+1)=Nh(k+1)+nhk
gamma(k+1,2)=gamma(k+1,2)+gk
else
Nh(k+1)=nhk
gamma(k+1,2)=gk
end if
end if
end do
!backward
do k=nbins,2,-1
nhk=Nh(k)
gk=gamma(k,2)
if ((nhk>0) .and. (nhk<30)) then
Nh(k)=-9
gamma(k,2)=0_dp
if (Nh(k-1)>0) then
Nh(k-1)=Nh(k-1)+nhk
gamma(k-1,2)=gamma(k-1,2)+gk
else
Nh(k-1)=nhk
gamma(k-1,2)=gk
end if
end if
end do
!
ni=0
! estimate semi-variance gamma(i,1)=h, gamma(i,2)=g(h)
do k=1,nbins
if (Nh(k) > 0) then
ni=ni+1
! Classsical
gamma(k,2)=gamma(k,2)/2.0_dp/real(Nh(k), dp)
!
! Cressi-Hawkins: adjust bias
!gamma(k,2)=0.5_dp/(0.457_dp+0.494_dp/dfloat(Nh(k)))*(gamma(k,2)/dfloat(Nh(k)))**4
!
gamma(k,1)=real(k, dp)*dh-real(ni, dp)*dh/2.0_dp
ni=0
else
ni=ni+1
gamma(k,1)=-9.0_dp
end if
end do
!
! scaling
if (flagMax) gmax=maxval(gamma, dim=1)
do k=1,nbins
if (gamma(k,1) > 0.0_dp) then
gamma(k,1) = gamma(k,1)/gmax(1)
gamma(k,2) = gamma(k,2)/gmax(2)
end if
end do
!
!keep variogram
end if
end subroutine EmpVar
end module mo_EmpVar
......@@ -103,12 +103,6 @@ contains
ii=Nk(i)
cell(k)%z(jd) = cell(k)%z(jd) + X(i) * MetSta(ii)%z(jd)
end do
! ! correct in case of negative values
! ! sick system of equations caused by lack of
! ! precipitation values in nearby stations
! if (correctNeg) then
! cell(k)%z = 0.0_dp
! end if
!
! only one station available /= 0
else if (ll /= nNmax .AND. nNmax == 1) then
......@@ -123,7 +117,6 @@ contains
! if all stations have the value 0
else if (ll == nNmax ) then
cell(k)%z(jd) = 0.0_dp
! avoid numerical instabilities --> IWD insverse weighted squared distance
! matrix of solver may become instable if only two or three stations are available
else if (ll /= nNmax .AND. nNmax == 2) then
......@@ -220,16 +213,17 @@ contains
! NOTE: only the upper triangular matrix is needed!
! call D_LSASF (A, B, X)
! print *, '<<<<<<<<<<<<<<'
! print *, 'rhs = ', B
! print *, 'rhs = ', B(:10)
! C = A
! print *, 'A = ', C(:, 1)
! X = B
! print *, 'A = ', C(:10, 1)
! print *, '<<<<<<<<<<<<<<'
call dgesv(size(A, 1), 1, A, size(A, 1), ipvt, B, size(A, 1), info)
! print *, '<<<<<<<<<<<<<<'
! print *, 'B = ', B
! print *, 'A = ', A(:, 1)
! print *, 'B = ', B(:10)
! print *, 'A = ', A(:10, 1)
! print *, '<<<<<<<<<<<<<<'
! print *, 'result: ', sum(C(:, 1) * B)
if (maxval(abs(matmul(C, B) - X)) .gt. 1e-14) print *, 'result: ', maxval(abs(matmul(C, B) - X))
X = B
if (info .ne. 0_i4) then
print *, '***WARNING: calculation of weights failed'
......@@ -238,9 +232,11 @@ contains
print *, '***WARNING: sum of weights is not 1, calculation of weights failed'
print *, 'sum of weights: ', sum(X(:nNmax))
end if
! print *, 'easting: ', cell%x
! print *, 'northing: ', cell%y
! print *, 'number of neighbors: ', nNmax
! print *, ''
! print *, 'ipvt: ', ipvt
! ! print *, 'ipvt: ', ipvt
! print *, 'info: ', info
! stop 'testing'
deallocate (A,B,ipvt)
......
......@@ -13,7 +13,7 @@ CONTAINS
use mo_kind, only: i4, sp, dp
use mo_netcdf, only: NcDataset, NcDimension, NcVariable
use mo_string_utils, only: num2str
use mainVar, only: gridMeteo, yStart
use mainVar, only: gridMeteo, yStart, mStart, dStart
use NetCDFVar, only: fileOut, author_name, variable_name
implicit none
......@@ -37,7 +37,9 @@ CONTAINS
! create variables
var_time = nc%setVariable('time', "i32", (/dim_time/))
! add some variable attributes
call var_time%setAttribute("units", "days since " // trim(num2str(yStart - 1, form='(I4)')) // "-12-31 12:00:00")
call var_time%setAttribute("units", "days since " // trim(num2str(yStart, form='(I4)')) // "-"// &
trim(num2str(mStart, form='(I0.2)')) // "-" // &
trim(num2str(dStart, form='(I0.2)')) // "-" // "12:00:00")
allocate(dummy(gridMeteo%ncols, gridMeteo%nrows))
do i = 1, gridMeteo%nrows
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment