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

Commit b85f3810 authored by Akash Koppa's avatar Akash Koppa
Browse files

added a feature to allow time buffering

parent 0880e492
File added
......@@ -40,7 +40,7 @@ yStart = 1990
mStart = 1
dStart = 1
yEnd = 1990
mEnd = 12
mEnd = 12
dEnd = 31
!
! -------------------------------------------------------------------
......
File added
ncols 6
nrows 10
xllcorner 4530000.0
yllcorner 5664000.0
cellsize 4000
NODATA_value -9999
......@@ -25,7 +25,7 @@ program ED_Kriging
use runControl , only: flagEDK, interMth, & ! flag for activate kriging, flag for 'OK' or 'EDK'
correctNeg, & ! pre or temp
flagVario ! flag for activate variogram estimation
use mainVar , only: yStart, yEnd, jStart, jEnd, & ! interpolation time periods
use mainVar , only: yStart, yEnd, jStart, jEnd, nSta, & ! interpolation time periods
grid, gridMeteo, & ! grid properties of input and output grid
nCell, MetSta, &
noDataValue
......@@ -52,13 +52,17 @@ program ED_Kriging
integer(i4) :: itimer
integer(i4) :: doy ! day of year
integer(i4) :: year, month, day ! current date
integer(i4) :: itime, itemp,sttemp,cnttemp ! loop variable for chunking write
integer(i4) :: jstarttmp, jendtmp ! more loop variables
!$ 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
integer(i4), allocatable :: Nk_old(:)
real(dp), allocatable :: X(:)
call print_start_message()
!$OMP PARALLEL
......@@ -107,6 +111,7 @@ program ED_Kriging
call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
call message('')
!write(*,*), "jStart = ",jStart
if (interMth .gt. 0) then
itimer = 4
......@@ -116,16 +121,72 @@ 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
! Akash----------------------------------------------------------------
do iCell = 1, nCell
! initialize cell
allocate(cell(iCell)%z(jStart:jEnd))
cell(iCell)%z = noDataValue
allocate(cell(iCell)%Nk_old(nSta))
cell(iCell)%Nk_old = -9999
end do
do iCell = 1, nCell
! initialize cell
allocate(cell(iCell)%W(nSta))
cell(iCell)%W = -9999
end do
if (mod((jEnd - jStart + 1),100) .eq. 0) then ! just use mod
iTime = ((jEnd - jStart + 1)/100)
else
iTime = ((jEnd - jStart + 1)/100) + 1
end if
write(*,*),"Total Number of Time Buffers = ",iTime
t = 0
bufferloop: do iTemp = 1, iTime
!call message('Time Loop Running')
!$OMP parallel default(shared) &
!$OMP private(iCell)
!$OMP do SCHEDULE(STATIC)
!call open_netcdf(nc_out, nc_data, nc_time) ! dont do this
!jStartTmp = jStart
!if (iTemp .lt. iTime) then
! jEndTmp = jStartTmp + 100
!else
! jEndTmp = (jEnd-jStart+1) - ((iTemp-1) * 100)
!end if ! use minimum to never exceed jEnd
jStartTmp = jStart + (iTemp - 1) * 100
if (iTemp .lt. iTime) then
jEndTmp = jStartTmp + 100 - 1
else
jEndTmp = jStartTmp + (jEnd-jStartTmp+1)
end if ! use minimum to never exceed jEnd
jEndTmp = min(jEndTmp, jEnd)
! write(*,*), "jStartTmp = ",jStartTmp," jEndTmp = ",jEndTmp,"jEnd = ",jEnd
do iCell = 1, nCell
! initialize cell ! deallocate similarly
allocate(cell(iCell)%z(jStartTmp:jEndTmp))
cell(iCell)%z = noDataValue
end do
!$OMP parallel default(shared) &
!$OMP private(iCell, X, Nk_old)
!$OMP do SCHEDULE(STATIC)
ncellsloop: do iCell = 1, nCell
! check DEM
......@@ -133,20 +194,46 @@ program ED_Kriging
cell(iCell)%z = gridMeteo%nodata_value
cycle
end if
!write(*,*),"Interpolation Method = ",interMth
!write(*,*),"Cell = ",iCell
!write(*,*),"jStartTmp = ",jStartTmp
!write(*,*),"jEndTmp = ",jEndTmp
!write(*,*),"dCS = ",dCS
!write(*,*),"MetSta = ",MetSta
!write(*,*),"Flag 1"
X = cell(iCell)%W
!write(*,*),"Flag 1.25"
!write(*,*),tempX
Nk_old = cell(iCell)%Nk_old
!write(*,*),"Flag 1.5"
!write(*,*),"X before call EDK: ",X
! interploation
select case (interMth)
case (1)
call EDK(iCell, jStart, jEnd, dCS, MetSta, dS, cell, doOK=.True.)
!call EDK(iCell, jStartTmp, jEndTmp, dCS, MetSta, dS, cell, doOK=.True.)
!call EDK(iCell, jStartTmp, jEndTmp, dCS, MetSta, dS, cell, tempX, tempNkOld, doOK=.True.)
call EDK(iCell, jStartTmp, jEndTmp, dCS, MetSta, dS, cell, X, Nk_old, doOK=.True.)
case (2)
call EDK(iCell, jStart, jEnd, dCS, MetSta, dS, cell)
!call EDK(iCell, jStartTmp, jEndTmp, dCS, MetSta, dS, cell)
!call EDK(iCell, jStartTmp, jEndTmp, dCS, MetSta, dS, cell, tempX, tempNkOld)
call EDK(iCell, jStartTmp, jEndTmp, dCS, MetSta, dS, cell, X, Nk_old)
end select
cell(iCell)%W = X
cell(iCell)%Nk_old = Nk_old
!write(*,*),"X after call EDK = ",X
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))
allocate(tmp_array(gridMeteo%ncols, gridMeteo%nrows, jEndTmp - jStartTmp + 1))
allocate(tmp_time(jEndTmp - jStartTmp + 1))
k = 0
do i = 1, gridMeteo%ncols
......@@ -167,15 +254,64 @@ program ED_Kriging
end do
end if
end do
t = 0
do i = 1, jEnd - jStart + 1
!t = 0
!do i = 1, jEnd - jStart + 1
! tmp_time(i) = t
! t = t + 1
!end do
do i = 1, jEndTmp - jStartTmp + 1
tmp_time(i) = t
t = t + 1
end do
call nc_time%setData(tmp_time)
call nc_data%setData(tmp_array)
deallocate(tmp_array, tmp_time)
!write(*,*),tmp_time
sttemp = nint(tmp_time(1)+1)
cnttemp = nint((tmp_time(size(tmp_time)) - sttemp))+2
!write(*,*),"array_size = ",size(tmp_array,1)," ",size(tmp_array,2)
!write(*,*),"array_size = ",size(tmp_array,3)
!write(*,*),"tmp_time_1 = ",sttemp
!write(*,*),"tmp_last = ",cnttemp
!write(*,*),"Start tmp_time",(/sttemp/)
!write(*,*),"End tmp_time",(/cnttemp/)
!write(*,*),"Start tmp_array",(/sttemp,1,1/)
!write(*,*),"End tmp_array",(/cnttemp,size(tmp_array,2),size(tmp_array,1)/)
!call nc_time%setData(values=tmp_time,start=(/sttemp/),cnt=(/cnttemp/))
!call nc_data%setData(values=tmp_array,start=(/sttemp,1,1/),cnt=(/cnttemp,size(tmp_array,2),size(tmp_array,1)/))
call nc_time%setData(values=tmp_time,start=(/sttemp/),cnt=(/cnttemp/))
call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
deallocate(tmp_array, tmp_time)
!deallocate(cell)
do iCell = 1, nCell
! initialize cell
deallocate(cell(iCell)%z)
!cell(iCell)%z = noDataValue
end do
! close netcdf if necessary
!call nc_out%close() ! outside
end do bufferloop
! Akash end----------------------------------------------------
do iCell = 1, nCell
! initialize cell
deallocate(cell(iCell)%Nk_old)
!cell(iCell)%z = noDataValue
end do
do iCell = 1, nCell
! initialize cell
deallocate(cell(iCell)%W)
!cell(iCell)%z = noDataValue
end do
! close netcdf if necessary
call nc_out%close()
......
......@@ -20,7 +20,7 @@ contains
! 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)
subroutine EDK(k, jStart, jEnd, dCS, MetSta, dS, cell, X, Nk_old, doOK)
use mo_kind, only : i4, dp, sp
use mainVar, only : MeteoStation, noDataValue, nSta, thresholdDist
use kriging, only : CellCoarser, dtoS
......@@ -29,25 +29,29 @@ contains
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
integer(i4), intent(in) :: k ! cell id
integer(i4), intent(in) :: jStart, jEnd
integer(i4), intent(inout) :: Nk_old(nSta) ! added Nk_old(nSta)
real(dp), intent(in) :: dCS(:, :) ! distances matrix
real(dp), allocatable, intent(inout) :: X(:) ! added X(:)
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(:,:)
integer(i4) :: Nk(nSta),nNmax !Nk_old !Nk(nSta), Nk_old(nSta), nNmax ! deleted: Nk_old(nSta)
real(dp), allocatable :: A (:,:), B(:), C(:,:) ! deleted: X(:)
real(dp), allocatable :: lamda(:)
real(dp) :: sumLamda
!write(*,*),"Flag 3"
!
! Check which stations have valid data ( =/ -9, 999 etc. ) store them in Nk(:)
! Check stations have valid data ( =/ -9, 999 etc. ) store them in Nk(:)
l = 0
ll = 0
Nk = 0
......@@ -56,7 +60,7 @@ contains
! switch ordinary kriging off if not explicitly given
doOK_loc = .False.
if (present(doOK)) doOK_loc = doOK
!write(*,*),"Flag 4"
! IF NK changed -> re-estimate weights
timeloop: do jd = jStart, jEnd
......@@ -81,7 +85,9 @@ contains
else
calc_weights = .True.
end if
!write(*,*),"Nk = ", Nk
!write(*,*),"Nk_old = ",Nk_old
!write(*,*),"calc_weights = ",calc_weights
!>>>> no value ! avoid indetermination
! avoid 0 value calculations ll == nNmax
......@@ -91,10 +97,12 @@ contains
if (calc_weights) then
! print *, 'cell: ', k
! print *, '#neighbors: ', nNmax
! print *, 'Nk: ', Nk
! print *, 'Nk_old: ', Nk_old
!print *, '#neighbors: ', nNmax
!print *, 'Nk: ', Nk
!print *, 'Nk_old: ', Nk_old
!write(*,*),"calc_weights inside = ",calc_weights
call get_kriging_weights(X, nNmax, Nk, doOK_loc, dCS(k, :), dS, cell(k), MetSta)
!write(*,*),"X after kriging weights: ",X
end if
! The BLUE of z is then:
......@@ -136,7 +144,7 @@ contains
end if
! stop 'TESTING'
end do timeloop
if (allocated(X)) deallocate (X)
!if (allocated(X)) deallocate (X)
!
! correct negative
if (correctNeg) cell(k)%z = merge(0._sp, cell(k)%z, (cell(k)%z .gt. -9999._sp) .and. (cell(k)%z .lt. 0.))
......@@ -155,7 +163,7 @@ contains
real(dp), allocatable, intent(out) :: X(:)
integer(i4), intent(in) :: nNmax
integer(i4), intent(in) :: Nk(nSta)
integer(i4), intent(in) :: Nk(nSta)!Nk(nSta)
logical, intent(in) :: doOK_loc
real(dp), intent(in) :: dCS(:) ! distances matrix
type(dtoS), intent(in) :: dS(:) ! distance among stations
......
......@@ -124,14 +124,17 @@ end module runControl
module kriging
use mo_kind, only: i4, sp, dp
use mainVar, only: nSta
real(dp) :: maxDist ! max distance [m] search stations
type CellCoarser
integer(i4) :: nNS ! No. Nearest Stations (NS) d <= maxDist
integer(i4), allocatable :: Nk_old(:) ! old stations (added)
integer(i4), dimension(:), allocatable :: listNS ! list of NS
real(dp) :: x ! x- coordinate
real(dp) :: y ! y- coordinate
real(dp) :: h ! (estimated) elevation [m] (from the nearest cells DEM)
real(sp), allocatable :: z(:) ! z values to be interpolated (OUTPUT)
real(dp), allocatable :: W(:)
end type CellCoarser
type(CellCoarser), dimension(:), allocatable :: cell ! EDK output
......@@ -141,8 +144,9 @@ module kriging
end type dtoS
type (dtoS), dimension(:), allocatable :: dS ! distance from Station i to all js
type (dtoS), dimension(:), allocatable :: dz2S ! squared diference of Z values
real(dp) :: xl, xr, yd, yu ! coordinates of the interpolation block
!
real(dp) :: xl, xr, yd, yu ! coordinates of the interpolation block
!real(dp) , allocatable :: X(:)
!
end module kriging
module VarFit
......
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