Feel free to join the next Helmholtz Hacky Hour #26 on Wednesday, April 21, 2021 from 2PM to 3PM!

Commit ed307296 authored by Friedrich Boeing's avatar Friedrich Boeing

Merge branch 'mem_optimization' into 'develop'

Memory optimization

See merge request !5
parents 255f1cc8 c92bef41
......@@ -3,11 +3,11 @@
# lines ahead concerning module independend builds
cmake_minimum_required(VERSION 3.5)
project(edk Fortran)
project(edk Fortran C CXX)
# The variable "CMAKE_BUILD_MODULE_SYSTEM_INDEPENDEND" can be set before executing cmake via a cache command:
# $cmake -DCMAKE_BUILD_MODULE_SYSTEM_INDEPENDEND:STRING=ON ..
# or cache file:
# or cache file:
# $cmake -C ../CMakeCacheFiles/eve ..
# or after executing CMake editing the CMakeCache.txt, preferably with a corresponding cmake editor i.e ccmake
set(CMAKE_BUILD_MODULE_SYSTEM_INDEPENDEND OFF CACHE STRING "build the module independend of the module system, so the build in the build tree works even after a module purge")
......
......@@ -14,7 +14,7 @@
! UPDATES
! Created Sa 21.03.2006
! Last Update Sa 11.06.2010 ! blocks, whole Germany
! Last Update Zi 04.02.2012 ! changed to general edk version
! Last Update Zi 04.02.2012 ! changed to general edk version
! ! (excluded block seperation)
!****************************************************************************
program ED_Kriging
......@@ -23,13 +23,13 @@ program ED_Kriging
use mo_print_message , only: print_start_message, print_end_message
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, tBuffer, nSta, DEMNcFlag, & ! interpolation time periods
grid, gridMeteo, & ! grid properties of input and output grid
nCell, MetSta, &
noDataValue
use kriging , only: dCS, dS, cell
grid, gridMeteo, & ! grid properties of input and output grid
nCell, MetSta, &
noDataValue
use kriging , only: edk_dist, cell
use mo_setVario , only: setVario, dMatrix
use mo_netcdf , only: NcDataset, NcVariable
use mo_write , only: open_netcdf
......@@ -38,7 +38,7 @@ program ED_Kriging
use mo_ReadData , only: readData
use NetCDFVar , only: invert_y
USE mo_timer, ONLY : &
timers_init, timer_start, timer_stop, timer_get ! Timing of processes
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
......@@ -63,9 +63,9 @@ program ED_Kriging
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(:)
! integer(i4), allocatable :: Nk_old(:)
! real(dp), allocatable :: X(:)
call print_start_message()
loop_factor = 10 ! factor for setting openMP loop size
......@@ -83,7 +83,7 @@ program ED_Kriging
itimer = 1
call timer_start(itimer)
call message('')
call message('')
call message(' >>> Reading data')
call message('')
call ReadData
......@@ -96,11 +96,11 @@ program ED_Kriging
!print *, 'nCell: ', nCell
!print *, "ncell_thread: ", ncell_thread
!print *, 'n_threads: ', n_threads
! print *, 'DEMNcFlag', DEMNcFlag
! print *, 'DEMNcFlag', DEMNcFlag
itimer = 2
call timer_start(itimer)
call message(' >>> Calculating distance matrix')
call message(' >>> Calculating distances and neighborhoods')
call message('')
! call distance matrix
call dMatrix
......@@ -109,21 +109,20 @@ program ED_Kriging
call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
call message('')
itimer = 3
call timer_start(itimer)
call message(' >>> Estimate variogram')
call message('')
! estimate variogram
call setVario(param)
! write variogram
! write variogram
if (flagVario) call WriteDataMeteo(0,0,2)
call timer_stop(itimer)
call message('')
call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
call message('')
!write(*,*), "jStart = ",jStart
!write(*,*), "jStart = ",jStart
if (interMth .gt. 0) then
itimer = 4
......@@ -133,165 +132,155 @@ 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)%Nk_old(nSta))
! cell(iCell)%Nk_old = nint(noDataValue)
! end do
if (mod((jEnd - jStart + 1),tBuffer) .eq. 0) then ! just use mod
iTime = ((jEnd - jStart + 1)/tBuffer)
else
iTime = ((jEnd - jStart + 1)/tBuffer) + 1
end if
write(*,*) "Total Number of Threads = ",n_threads
write(*,*) "Total Number of Time Buffers = ",iTime
t = 0
bufferloop: do iTemp = 1, iTime
write(*,*) " >>> Started buffer #", iTemp
jStartTmp = jStart + (iTemp - 1) * tBuffer
if (iTemp .lt. iTime) then
jEndTmp = jStartTmp + tBuffer - 1
else
jEndTmp = jStartTmp + (jEnd-jStartTmp+1)
end if ! use minimum to never exceed jEnd
jEndTmp = min(jEndTmp, jEnd)
do iCell = 1, nCell
! initialize cell ! deallocate similarly
allocate(cell(iCell)%z(jStartTmp:jEndTmp))
cell(iCell)%z = real(noDataValue, sp)
end do
!print *, iTemp, iTime
!$OMP parallel default(shared) &
!$OMP private(iThread, iCell)
!$OMP do SCHEDULE(dynamic)
do iThread = 1, loop_factor * n_threads
! print *, 'thread: ', iThread, " start"
ncellsloop: do iCell = (iThread - 1) * ncell_thread + 1, min(iThread * ncell_thread, ncell)
! check DEM
if (nint(cell(iCell)%h) == grid%nodata_value ) then
cell(iCell)%z = gridMeteo%nodata_value
cycle
end if
! interploation
call EDK(iCell, jStartTmp, jEndTmp, edk_dist, MetSta, cell(iCell), doOK=(interMth==1))
end do ncellsloop
! print *, 'thread: ', iThread, " end"
end do
!$OMP end do
!$OMP end parallel
if (DEMNcFlag == 1) then
! write output
allocate(tmp_array(gridMeteo%nrows, gridMeteo%ncols, jEndTmp - jStartTmp + 1))
allocate(tmp_time(jEndTmp - jStartTmp + 1))
k = 0
if (invert_y) then
do i = 1, gridMeteo%ncols
do j = 1, gridMeteo%nrows
k = k + 1
tmp_array(j,gridMeteo%ncols - i + 1,:) = cell(k)%z
end do
end do
else
do i = 1, gridMeteo%ncols
do j = 1, gridMeteo%nrows
k = k + 1
tmp_array(j,i,:) = cell(k)%z
end do
end do
end if
do iCell = 1, nCell
! initialize cell
allocate(cell(iCell)%Nk_old(nSta))
cell(iCell)%Nk_old = nint(noDataValue)
end do
if (mod((jEnd - jStart + 1),tBuffer) .eq. 0) then ! just use mod
iTime = ((jEnd - jStart + 1)/tBuffer)
else
iTime = ((jEnd - jStart + 1)/tBuffer) + 1
end if
write(*,*),"Total Number of Time Buffers = ",iTime
t = 0
bufferloop: do iTemp = 1, iTime
jStartTmp = jStart + (iTemp - 1) * tBuffer
if (iTemp .lt. iTime) then
jEndTmp = jStartTmp + tBuffer - 1
else
jEndTmp = jStartTmp + (jEnd-jStartTmp+1)
end if ! use minimum to never exceed jEnd
jEndTmp = min(jEndTmp, jEnd)
do iCell = 1, nCell
! initialize cell ! deallocate similarly
allocate(cell(iCell)%z(jStartTmp:jEndTmp))
cell(iCell)%z = noDataValue
do i = 1, jEndTmp - jStartTmp + 1
tmp_time(i) = t
t = t + 1
end do
!print *, iTemp, iTime
!$OMP parallel default(shared) &
!$OMP private(iThread, iCell, X, Nk_old)
!$OMP do SCHEDULE(dynamic)
do iThread = 1, loop_factor * n_threads
! print *, 'thread: ', iThread, " start"
ncellsloop: do iCell = (iThread - 1) * ncell_thread + 1, min(iThread * ncell_thread, 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 EDK(iCell, jStartTmp, jEndTmp, dCS, MetSta, dS, cell, cell(iCell)%W, cell(iCell)%Nk_old, doOK=.True.)
case (2)
call EDK(iCell, jStartTmp, jEndTmp, dCS, MetSta, dS, cell, cell(iCell)%W, cell(iCell)%Nk_old)
end select
end do ncellsloop
! print *, 'thread: ', iThread, " end"
sttemp = nint(tmp_time(1)+1)
cnttemp = nint((tmp_time(size(tmp_time)) - sttemp))+2
!write(*,*),"Final Output ",shape(tmp_array)
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/))
call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
else
! write output
allocate(tmp_array(gridMeteo%ncols, gridMeteo%nrows, jEndTmp - jStartTmp + 1))
allocate(tmp_time(jEndTmp - jStartTmp + 1))
k = 0
do i = 1, gridMeteo%ncols
! do j = 1, gridMeteo%nrows
! k = k + 1
! tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
! end do
! end do
if (invert_y) then
do j = gridMeteo%nrows, 1, -1
k = k + 1
tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
end do
else
do j = 1, gridMeteo%nrows
k = k + 1
tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
end do
end if
end do
!$OMP end do
!$OMP end parallel
if (DEMNcFlag == 1) then
! write output
allocate(tmp_array(gridMeteo%nrows, gridMeteo%ncols, jEndTmp - jStartTmp + 1))
allocate(tmp_time(jEndTmp - jStartTmp + 1))
k = 0
if (invert_y) then
do i = 1, gridMeteo%ncols
do j = 1, gridMeteo%nrows
k = k + 1
tmp_array(j,gridMeteo%ncols - i + 1,:) = cell(k)%z
end do
end do
else
do i = 1, gridMeteo%ncols
do j = 1, gridMeteo%nrows
k = k + 1
tmp_array(j,i,:) = cell(k)%z
end do
end do
end if
do i = 1, jEndTmp - jStartTmp + 1
tmp_time(i) = t
t = t + 1
end do
sttemp = nint(tmp_time(1)+1)
cnttemp = nint((tmp_time(size(tmp_time)) - sttemp))+2
!write(*,*),"Final Output ",shape(tmp_array)
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/))
call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
else
! write output
allocate(tmp_array(gridMeteo%ncols, gridMeteo%nrows, jEndTmp - jStartTmp + 1))
allocate(tmp_time(jEndTmp - jStartTmp + 1))
k = 0
do i = 1, gridMeteo%ncols
! do j = 1, gridMeteo%nrows
! k = k + 1
! tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
! end do
! end do
if (invert_y) then
do j = gridMeteo%nrows, 1, -1
k = k + 1
tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
end do
else
do j = 1, gridMeteo%nrows
k = k + 1
tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
end do
end if
end do
!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
!write(*,*),tmp_time
sttemp = nint(tmp_time(1)+1)
cnttemp = nint((tmp_time(size(tmp_time)) - sttemp))+2
!write(*,*),"Final Output ",shape(tmp_array)
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/))
call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
end if
!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
!write(*,*),tmp_time
sttemp = nint(tmp_time(1)+1)
cnttemp = nint((tmp_time(size(tmp_time)) - sttemp))+2
!write(*,*),"Final Output ",shape(tmp_array)
deallocate(tmp_array, tmp_time)
!deallocate(cell)
do iCell = 1, nCell
! initialize cell
deallocate(cell(iCell)%z)
!cell(iCell)%z = noDataValue
end do
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/))
call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
end if
deallocate(tmp_array, tmp_time)
!deallocate(cell)
! close netcdf if necessary
!call nc_out%close() ! outside
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
end do bufferloop
! close netcdf if necessary
call nc_out%close()
......@@ -307,4 +296,3 @@ program ED_Kriging
call print_end_message()
!
end program ED_Kriging
!*************************************************************************
!EmpVar.f90
!EmpVar.f90
!AUTHOR: Luis Samaniego-eguiguren
!
!
!DESCRIPTION
! This program calculates the empirical variogram
! This program calculates the empirical variogram
! gamma(:,1) : distance
! gamme(:,2) : semivarance
! botm are normalized with gmax(:)
......@@ -31,24 +31,7 @@ subroutine EmpVar(jd, flagMax)
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
......@@ -64,7 +47,7 @@ subroutine EmpVar(jd, flagMax)
!
! selecting pair - bins
if (jd == jStart) then
nbins = ceiling(hMax/dh)
nbins = ceiling(hMax/dh)
if (.not. allocated (nH)) allocate (Nh(nbins), gamma(nbins,2))
Nh = 0
gamma = 0.0_dp
......@@ -73,20 +56,20 @@ subroutine EmpVar(jd, flagMax)
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
if (edk_dist%getZ2(i,j,jd) /= noDataValue ) then
! take values up to max distance
if (dS(i)%S(j) > hMax ) cycle
if (edk_dist%getSS(i,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
! if (flagVarTyp == 2 ) then
if ( dabs( MetSta(i)%h - MetSta(j)%h ) / edk_dist%getSS(i,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))
k=max(1, ceiling(edk_dist%getSS(i,j)/dh))
Nh(k)=Nh(k)+1
gamma(k,2)=gamma(k,2) + dz2S(i)%S(j)
gamma(k,2)=gamma(k,2) + edk_dist%getZ2(i,j,jd)
end if
end do
end do
......@@ -102,7 +85,7 @@ subroutine EmpVar(jd, flagMax)
if ((nhk>0) .and. (nhk<30)) then
Nh(k)=-9
gamma(k,2)=0.0_dp
if (Nh(k+1)>0) then
if (Nh(k+1)>0) then
Nh(k+1)=Nh(k+1)+nhk
gamma(k+1,2)=gamma(k+1,2)+gk
else
......@@ -118,7 +101,7 @@ subroutine EmpVar(jd, flagMax)
if ((nhk>0) .and. (nhk<30)) then
Nh(k)=-9
gamma(k,2)=0_dp
if (Nh(k-1)>0) then
if (Nh(k-1)>0) then
Nh(k-1)=Nh(k-1)+nhk
gamma(k-1,2)=gamma(k-1,2)+gk
else
......@@ -138,7 +121,7 @@ subroutine EmpVar(jd, flagMax)
!
! 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
......@@ -147,9 +130,9 @@ subroutine EmpVar(jd, flagMax)
end if
end do
!
! scaling
! 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)
......
......@@ -14,7 +14,7 @@ module mo_ReadData
real(dp) :: right
real(dp) :: bottom
end type extend
contains
!*************************************************************************
! SUBROUTINE Read Database Precipitation
......@@ -37,7 +37,7 @@ contains
distZero, flagVario, fNameVario, flagEDK
use NetCDFVar, only: FileOut, author_name, projection_name, variable_name, variable_unit, &
variable_standard_name, variable_calendar_type, &
variable_long_name, ncIn_variable_name, ncIn_dem_variable_name, &
variable_long_name, ncIn_variable_name, ncIn_dem_variable_name, &
ncIn_yCoord_name, ncIn_xCoord_name, invert_y, &
ncOut_dem_variable_name, ncOut_dem_yCoord_name, ncOut_dem_xCoord_name, &
ncOut_dem_Latitude, ncOut_dem_Longitude
......@@ -61,7 +61,7 @@ contains
DataPathOut, FileOut, fNameSTA, cellFactor, DataConvertFactor, OffSet, InterMth, correctNeg, &
distZero, author_name, projection_name, variable_name, variable_unit, variable_long_name, &
variable_standard_name, variable_calendar_type, &
yStart, mStart, dStart, yEnd, mEnd, dEnd,tBuffer, maxDist, flagVario, vType, nParam, &
yStart, mStart, dStart, yEnd, mEnd, dEnd,tBuffer, maxDist, flagVario, vType, nParam, &
fNameVario, dh, hMax, ncIn_variable_name, ncIn_dem_variable_name, ncIn_yCoord_name, &
ncIn_xCoord_name, invert_y, &
ncOut_dem_variable_name, ncOut_dem_yCoord_name, ncOut_dem_xCoord_name, ncOut_dem_Latitude,&
......@@ -117,11 +117,11 @@ contains
print*, "***ERROR: Variogram type not known: ", vType
end if
write(*,*) 'variogram parameters: '
write(*,'(3(a12), a18)') 'nugget', 'sill', 'range', 'Variogramtype'
write(*,*) 'variogram parameters: '
write(*,'(3(a12), a18)') 'nugget', 'sill', 'range', 'Variogramtype'
write(*,'(3(f12.2), i18)') (beta(i), i=1,nParam), vType
! determine whether DEM is an nc file or not
! determine whether DEM is an nc file or not
call divide_string(fNameDEM, '.', fNameDEM_parts)
isDEMNcFile = (fNameDEM_parts(size(fNameDEM_parts,1)) .eq. 'nc')
deallocate(fNameDEM_parts)
......@@ -134,7 +134,7 @@ contains
else
call ReadDEM
end if
! determine whether input path is actually a nc file
call divide_string(DataPathIn, '.', DataPathIn_parts)
isNcFile = (DataPathIn_parts(size(DataPathIn_parts,1)) .eq. 'nc')
......@@ -147,7 +147,7 @@ contains
else
! read look up table for meteorological stations
call ReadStationLut
! read whole METEO data
call ReadDataMeteo
end if
......@@ -160,8 +160,8 @@ contains
subroutine ReadStationLut
use mo_kind, only: i4, dp
use mainVar, only: nSta, MetSta, grid