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

Commit 4e2f5306 authored by Stephan Thober's avatar Stephan Thober
Browse files

Merge branch 'interpolation' into 'develop'

Interpolation

See merge request !3
parents 2ecdb269 a6c59ff6
!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