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

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

started interpolation

- finished read of meteo data from netcdf file
- bug in calculation of weights
parent 2ecdb269
!EDK-Grid***Main configuration file for Precipitation for GER***
&mainVars ! namelist mainpaths
flagMthTyp = 1
flagVarTyp = 1
!
! Correct interpolated negative values to zero
correctNeg = .False.
!
! Set values further away than the distance threshold to zero
distZero = .True.
noDataValue = -9
!
! DataPathIn contains the path to the meteorological input files
! it can be a netcdf file, then a 3-dim field is expected and
! further coordinate and variable names are expected
DataPathIn = "check/pre_data/"
! DataPathIn = "output/input.nc"
!
! Specifications of variable and coordinate names in case DataPathIn
! is a netcdf file
ncIn_variable_name = "pre"
!
! name of coordinates, these must be in [m]
ncIn_yCoord_name = "northing"
ncIn_xCoord_name = "easting"
!
!
fNameDEM = "check/dem.asc"
DataPathOut = "output/"
FileOut = "test.nc"
FileOut = "edk.nc"
fNameSTA = "check/pre_data/Stations_in_study_domain.txt"
cellFactor = 40
DataConvertFactor = 1.0e-1
flagEDK = .TRUE.
!
! SET INTERPOLATION
!
! InterMth = 2 -> EDK
! InterMth = 1 -> OK
! InterMth = 0 -> No interpolation
InterMth = 2
!
! PROCESSING PERIOD
yStart = 1990
mStart = 1
dStart = 1
yEnd = 1990
mEnd = 12
dEnd = 31
!
! maximum search distance for interpolation [m]
maxDist = 1.5e5
!
! estimate variogram
flagVario = .False.
vType = 1
!
! number of variogram parameters
nParam = 3
!
! file name where to store the variogram parameters
fNameVario = "output/varFit.txt"
dh = 1.0e3 ! binsize for variogram [m]
hMax = 1.5e5 ! max distance h [m]
hMax = 1.5e5 ! max distance h for variogram estimation [m]
! nc output specification
author_name = 'Stephan Thober'
variable_name = 'pre'
......
......@@ -20,13 +20,14 @@
program ED_Kriging
use mo_kind , only: i4, dp
use mo_julian , only: NDAYS, NDYIN
use mo_julian , only: NDAYS, NDYIN, dec2date, julday
use runControl , only: flagEDK, interMth, & ! flag for activate kriging, flag for 'OK' or 'EDK'
correctNeg, & ! pre or temp
flagVario ! flag for activate variogram estimation
correctNeg, & ! pre or temp
flagVario ! flag for activate variogram estimation
use mainVar , only: yStart, yEnd, jStart, jEnd, & ! interpolation time periods
grid, gridMeteo, & ! grid properties of input and output grid
nCell, MetSta
grid, gridMeteo, & ! grid properties of input and output grid
nCell, MetSta, &
noDataValue
use kriging , only: dCS, dS
use mo_setVario , only: setVario, dMatrix
use mo_netcdf , only: NcDataset, NcVariable
......@@ -34,26 +35,29 @@ program ED_Kriging
use mo_message , only: message
use kriging
use mo_EDK , only: EDK
use mo_ReadData , only: readData
USE mo_timer, ONLY : &
timers_init, timer_start, timer_stop, timer_get ! Timing of processes
!$ use mo_string_utils, ONLY : num2str
timers_init, timer_start, timer_stop, timer_get ! Timing of processes
use mo_string_utils, ONLY : num2str
!$ use omp_lib, ONLY : OMP_GET_NUM_THREADS ! OpenMP routines
implicit none
character(256) :: fname
character(256) :: author_name
character(256) :: vname_data
integer(i4) :: icell ! loop varaible for cells
integer(i4) :: jDay ! loop variable - current julian day
integer(i4) :: itimer
integer(i4) :: doy ! day of year
integer(i4) :: year, month, day ! current date
!$ integer(i4) :: n_threads ! OpenMP number of parallel threads
real(dp), dimension(:,:), allocatable :: tmp_array ! temporal array for output
real(dp) :: param(3) ! variogram parameters
type(NcDataset) :: nc_out
type(NcVariable) :: nc_data, nc_time
character(256) :: fname
character(256) :: author_name
character(256) :: vname_data
integer(i4) :: i, j, k
integer(i4) :: icell ! loop varaible for cells
integer(i4) :: jDay ! loop variable - current julian day
integer(i4) :: itimer
integer(i4) :: doy ! day of year
integer(i4) :: year, month, day ! current date
!$ integer(i4) :: n_threads ! OpenMP number of parallel threads
real(sp), allocatable :: tmp_array(:, :, :) ! temporal array for output
real(sp), allocatable :: tmp_time(:) ! temporal array for time output
real(dp) :: param(3) ! variogram parameters
type(NcDataset) :: nc_out
type(NcVariable) :: nc_data, nc_time
!$OMP PARALLEL
!$ n_threads = OMP_GET_NUM_THREADS()
......@@ -70,16 +74,7 @@ program ED_Kriging
call message('')
call message(' >>> Reading data')
call message('')
call ReadDataMain
! read DEM
call ReadDEM
! read look up table for meteorological stations
call ReadStationLut
! read whole METEO data
call ReadDataMeteo
call ReadData
call timer_stop(itimer)
call message('')
call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
......@@ -111,7 +106,6 @@ program ED_Kriging
call message('')
if (interMth .gt. 0) then
itimer = 4
call timer_start(itimer)
......@@ -120,50 +114,49 @@ program ED_Kriging
! open netcdf if necessary
call open_netcdf(nc_out, nc_data, nc_time)
timeloop: do jday = jStart, jEnd
call NDYIN(jday, day, month, year)
doy = jday - NDAYS(1,1,year) + 1
print *, 'YEAR: ',year, 'DOY: ', doy
!$OMP parallel default(shared) &
!$OMP private(iCell)
!$OMP do SCHEDULE(STATIC)
ncellsloop: do iCell = 1, nCell
! check DEM
if (nint(cell(iCell)%h) == grid%nodata_value ) then
cell(iCell)%z = gridMeteo%nodata_value
cycle
end if
! interploation
select case (interMth)
case (1)
call OK(jday, iCell)
case (2)
call EDK(jday, iCell, dCS, MetSta, dS, cell)
end select
end do ncellsloop
!$OMP end do
!$OMP end parallel
! correct negative values
if (correctNeg) then
where ((cell(:)%z .LT. 0.0_sp) .AND. (cell(:)%z .GT. real(grid%nodata_value, sp)) )
cell(:)%z = 0.0_sp
end where
end if
! write output
allocate(tmp_array(gridMeteo%nrows, gridMeteo%ncols)); tmp_array=real(grid%nodata_value, dp)
tmp_array = real(reshape(cell(:)%z,(/gridMeteo%nrows, gridMeteo%ncols/)), dp)
call nc_time%setData(jday-jStart + 1, start=(/jday - jStart + 1/))
call nc_data%setData(tmp_array, start=(/1, 1, jday - jStart + 1/))
!$OMP parallel default(shared) &
!$OMP private(iCell)
!$OMP do SCHEDULE(STATIC)
ncellsloop: do iCell = 1, nCell
deallocate(tmp_array)
! initialize cell
allocate(cell(iCell)%z(jStart:jEnd))
cell(iCell)%z = noDataValue
end do timeloop
! check DEM
if (nint(cell(iCell)%h) == grid%nodata_value ) then
cell(iCell)%z = gridMeteo%nodata_value
cycle
end if
! interploation
select case (interMth)
case (1)
call EDK(iCell, jStart, jEnd, dCS, MetSta, dS, cell, doOK=.True.)
case (2)
call EDK(iCell, jStart, jEnd, dCS, MetSta, dS, cell)
end select
end do ncellsloop
!$OMP end do
!$OMP end parallel
! write output
allocate(tmp_array(gridMeteo%ncols, gridMeteo%nrows, jEnd - jStart + 1))
allocate(tmp_time(jEnd - jStart + 1))
k = 0
do i = 1, gridMeteo%ncols
do j = 1, gridMeteo%nrows
k = k + 1
tmp_array(i, j, :) = cell(k)%z
end do
end do
do i = 1, jEnd - jStart + 1
tmp_time(i) = i
end do
call nc_time%setData(tmp_time)
call nc_data%setData(tmp_array)
deallocate(tmp_array, tmp_time)
! close netcdf if necessary
call nc_out%close()
......@@ -175,9 +168,9 @@ program ED_Kriging
end if
! deallocate memory
call clean
! very important for check cases
write(*,*) 'Kriging finished!'
!
end program
end program ED_Kriging
!****************************************************************************
!
! SUBROUTINE: Ordinary Kriging
!
! PURPOSE: Perform OK (daily)
! Output: output gridsize = cellFactor * gridsize DEM
!
! Created L. Samaniego 14.04.2006
! Last Update Sa 14.04.2006
!
!****************************************************************************
subroutine OK(jd,k)
use mainVar
use kriging
use runControl
use varfit, only : beta
use mo_kind, only : i4, dp
use mo_setVario, only : tVar
! use LSASF_INT
! use lapack95, only: gesv
implicit none
integer(i4), intent(in) :: jd ! day
integer(i4), intent(in) :: k ! cell id
integer(i4) :: i, j, l, ll
integer(i4) :: ii, jj
integer(i4) :: Nk(nSta), nNmax
real(dp), allocatable :: A (:,:), B(:), X(:)
real(dp), allocatable :: lamda(:)
real(dp) :: sumLamda
!
! check DEM
if (cell(k)%h == grid%nodata_value ) then
cell(k)%z = gridMeteo%noData_value
return
end if
!
! Check which stations have valid data ( =/ -9, 999 etc. ) store them in Nk(:)
l = 0
do i= 1, cell(k)%nNS
j = cell(k)%listNS(i)
if ( MetSta(j)%z(jd) /= noDataValue ) then
l = l + 1
Nk(l) = j
end if
end do
nNmax = l
!
!>>>> no value ! avoid indetermination
! avoid 0 value calculations ll == nNmax
! avoid calculations where only 1 station is available
! avoid numerical instabilities nNmax == 2 (may happen that the solver matrix becomes singular)
if (.not. ( ll == nNmax .or. nNmax == 1 .or. nNmax == 2 ) ) then
! initialize matrices
allocate (A(nNmax+1,nNmax+1), B(nNmax+1), X(nNmax+1))
!
! assembling the system of equations: OK
A=0.0_dp
do i=1,nNmax-1
ii=Nk(i)
do j=i+1,nNmax
jj=Nk(j)
! available only the upper triangular matrix
if (jj > ii) then
A(i,j)=tVar(dS(ii)%S(jj),beta(1),beta(2),beta(3))
else
A(i,j)=tVar(dS(jj)%S(ii),beta(1),beta(2),beta(3))
end if
end do
A(i,nNmax+1) = 1.0_dp
B(i)=tVar(dCS(k,ii),beta(1),beta(2),beta(3))
end do
A(nNmax,nNmax+1) = 1.0_dp
ii=Nk(nNmax)
B(nNmax)=tVar(dCS(k,ii),beta(1),beta(2),beta(3))
!
B(nNmax+1) = 1.0_dp
!NOTE: only the upper triangular matrix is needed!
! call D_LSASF (A, B, X)
call dgesv(A, B)
X = B
!
! The BLUE of z is then:
cell(k)%z = 0.
do i=1,nNmax
ii=Nk(i)
cell(k)%z = cell(k)%z + real(X(i) * MetSta(ii)%z(jd))
end do
!
! only one station available /= 0
else if (ll /= nNmax .AND. nNmax == 1) then
ii = Nk(1)
cell(k)%z = MetSta(ii)%z(jd)
!
! for precipitation, distant values are set to zero
if (distZero .and. dCS(k,ii) .gt. thresholdDist) then
cell(k)%z = 0.0_dp
end if
!
! if all stations have the value 0
! if all stations have the value 0
else if (ll == nNmax ) then
cell(k)%z = 0.0_dp
! avoid numerical instabilities --> IWD insverse weighted squared distance
! matrix of DLSASF may become instable if only two or three stations are available
else if (ll /= nNmax .AND. nNmax == 2) then
allocate (lamda(nNmax))
do i=1,nNmax
ii=Nk(i)
lamda(i)=1/dCS(k,ii)/dCS(k,ii)
end do
sumLamda=sum(lamda)
lamda=lamda/sum(lamda)
cell(k)%z = 0.0_dp
do i=1,nNmax
ii=Nk(i)
cell(k)%z = cell(k)%z + lamda(i) * MetSta(ii)%z(jd)
end do
deallocate (lamda)
end if
end subroutine OK
This diff is collapsed.
......@@ -7,70 +7,180 @@ module mo_EDK
public :: EDK
contains
!****************************************************************************
!
! SUBROUTINE: External Drift Kriging
!
! PURPOSE: Perform EDK (daily)
! Output: output gridsize = cellFactor * gridsize DEM
! UPDATES
! Created L. Samaniego 22.03.2006
! Last Update Sa 14.04.2006
! Last Update Sa 14.07.2010 DEM ckeck
! Last Update Zi 05.05.2011 IWD if No stations < 2
! Last Update Zi 07.02.2012 correct prec data < 0 --> = 0
!****************************************************************************
subroutine EDK(jd, k, dCS, MetSta, dS, cell)
use mo_kind, only : i4, dp
use mainVar, only : MeteoStation, noDataValue, nSta, thresholdDist
use kriging, only : CellCoarser, dtoS
use runControl, only : correctNeg, distZero
use varfit, only : beta
use mo_setVario, only : tVar
implicit none
! input / output variables
integer(i4), intent(in) :: jd ! julian day
integer(i4), intent(in) :: k ! cell id
real(dp), intent(in) :: dCS(:, :) ! distances matrix
type(MeteoStation), intent(in) :: MetSta(:) ! MeteoStation input
type(dtoS), intent(in) :: dS(:) ! distance among stations
type(CellCoarser), intent(inout) :: cell(:) ! cell specification
! local variables
integer(i4) :: i, j, l, ll, info
integer(i4) :: ii, jj
integer(i4) :: Nk(nSta), nNmax
integer(i4), allocatable :: ipvt(:)
real(dp), allocatable :: A (:,:), B(:), X(:) !, C(:,:)
real(dp), allocatable :: lamda(:)
real(dp) :: sumLamda
!****************************************************************************
!
! Check which stations have valid data ( =/ -9, 999 etc. ) store them in Nk(:)
l = 0
ll = 0
do i = 1, cell(k)%nNS
j = cell(k)%listNS(i)
if ( MetSta(j)%z(jd) /= noDataValue ) then
!*** zero field ***********************************
if (MetSta(j)%z(jd) == 0.0_dp ) ll = ll + 1
!**********************************************************
l = l + 1
Nk(l) = j
end if
end do
nNmax = l
!>>>> no value ! avoid indetermination
! avoid 0 value calculations ll == nNmax
! avoid calculations where only 1 station is available
! avoid numerical instabilities nNmax == 2 (may happen that the solver matrix becomes singular)
if (.not. ( ll == nNmax .or. nNmax == 1 .or. nNmax == 2 ) ) then
! SUBROUTINE: External Drift Kriging
!
! PURPOSE: Perform EDK (daily)
! Output: output gridsize = cellFactor * gridsize DEM
! UPDATES
! Created L. Samaniego 22.03.2006
! Last Update Sa 14.04.2006
! Last Update Sa 14.07.2010 DEM ckeck
! Last Update Zi 05.05.2011 IWD if No stations < 2
! Last Update Zi 07.02.2012 correct prec data < 0 --> = 0
!****************************************************************************
subroutine EDK(k, jStart, jEnd, dCS, MetSta, dS, cell, doOK)
use mo_kind, only : i4, dp, sp
use mainVar, only : MeteoStation, noDataValue, nSta, thresholdDist
use kriging, only : CellCoarser, dtoS
use runControl, only : correctNeg, distZero
implicit none
! input / output variables
integer(i4), intent(in) :: k ! cell id
integer(i4), intent(in) :: jStart, jEnd
real(dp), intent(in) :: dCS(:, :) ! distances matrix
type(MeteoStation), intent(in) :: MetSta(:) ! MeteoStation input
type(dtoS), intent(in) :: dS(:) ! distance among stations
type(CellCoarser), intent(inout) :: cell(:) ! cell specification
logical, optional, intent(in) :: doOK ! switch do ordinary kriging
! local variables
logical :: doOK_loc, calc_weights
integer(i4) :: jd ! julian day
integer(i4) :: i, j, l, ll
integer(i4) :: ii, jj
integer(i4) :: Nk(nSta), Nk_old(nSta), nNmax
real(dp), allocatable :: A (:,:), B(:), X(:) , C(:,:)
real(dp), allocatable :: lamda(:)
real(dp) :: sumLamda
!
! Check which stations have valid data ( =/ -9, 999 etc. ) store them in Nk(:)
l = 0
ll = 0
Nk = 0
Nk_old = 0
! switch ordinary kriging off if not explicitly given
doOK_loc = .False.
if (present(doOK)) doOK_loc = doOK
! IF NK changed -> re-estimate weights
timeloop: do jd = jStart, jEnd
Nk_old = Nk
Nk = 0_i4
l = 0
ll = 0
do i = 1, cell(k)%nNS
j = cell(k)%listNS(i)
if ( MetSta(j)%z(jd) /= noDataValue ) then
!*** zero field ***********************************
if (MetSta(j)%z(jd) == 0.0_dp ) ll = ll + 1
!**********************************************************
l = l + 1
Nk(l) = j
end if
end do
nNmax = l
if (all(Nk == Nk_old) .and. allocated(X)) then
calc_weights = .False.
else
calc_weights = .True.
end if
!>>>> no value ! avoid indetermination
! avoid 0 value calculations ll == nNmax
! avoid calculations where only 1 station is available
! avoid numerical instabilities nNmax == 2 (may happen that the solver matrix becomes singular)
if (.not. ( ll == nNmax .or. nNmax == 1 .or. nNmax == 2 ) ) then
if (calc_weights) then
! print *, 'cell: ', k
! print *, '#neighbors: ', nNmax
! print *, 'Nk: ', Nk
! print *, 'Nk_old: ', Nk_old
call get_kriging_weights(X, nNmax, Nk, doOK_loc, dCS(k, :), dS, cell(k), MetSta)
end if
! The BLUE of z is then:
cell(k)%z(jd) = 0.
do i=1,nNmax
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
ii = Nk(1)
cell(k)%z(jd) = MetSta(ii)%z(jd)
!
! for precipitation, distant values are set to zero
if (distZero .and. dCS(k,ii) .gt. thresholdDist) then
cell(k)%z(jd) = 0.0_dp
end if
!
! 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
allocate (lamda(nNmax))
do i=1,nNmax
ii=Nk(i)
lamda(i)=1/dCS(k,ii)/dCS(k,ii)
end do
sumLamda=sum(lamda)
lamda=lamda/sum(lamda)
cell(k)%z = 0.0_dp
do i=1,nNmax
ii=Nk(i)
cell(k)%z(jd) = cell(k)%z(jd) + lamda(i) * MetSta(ii)%z(jd)
end do
deallocate (lamda)
end if
! stop 'TESTING'
end do timeloop
if (allocated(X)) deallocate (X)
!
! correct negative