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

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

Merge branch 'feature_openMP' into 'develop'

Feature open mp

See merge request !2
parents 6a932321 b863aed8
......@@ -96,7 +96,7 @@ SHELL = /bin/bash
# . is current directory, .. is parent directory
SRCPATH := ./src ./lib # where are the source files; use test_??? to run a test directory
PROGPATH := . # where shall be the executable
CONFIGPATH := /home/thober/lib/makefile_chs/make.config # where are the $(system).$(compiler) files
CONFIGPATH := /home/thober/lib/mhm/develop/make.config # where are the $(system).$(compiler) files
MAKEDPATH := $(CONFIGPATH) # where is the make.d.sh script
CHECKPATH := . # path for $(CHECKPATH)/test* and $(CHECKPATH)/check* directories if target is check
DOXCONFIG := ./doxygen.config # the doxygen config file
......@@ -106,24 +106,24 @@ LIBNAME := #libminpack.a # Name of library
#
# Options
# Systems: eve and personal computers such as mcimac for Matthias Cuntz' iMac; look in $(MAKEDPATH) or type 'make info'
system := eve
system := stmac
# Compiler: intelX, gnuX, nagX, sunX, where X stands for version number, e.g. intel13;
# look at $(MAKEDPATH)/$(system).alias for shortcuts or type 'make info'
compiler := intel
compiler := gnu71
# Releases: debug, release
release := release
release := debug
# Netcdf versions (Network Common Data Form): netcdf3, netcdf4, [anything else]
netcdf := netcdf4
# LAPACK (Linear Algebra Pack): true, [anything else]
lapack :=
lapack := true
# MKL (Intel's Math Kernel Library): mkl, mkl95, [anything else]
mkl := mkl95
mkl :=
# Proj4 (Cartographic Projections Library): true, [anything else]
proj :=
# IMSL (IMSL Numerical Libraries): vendor, imsl, [anything else]
imsl :=
# OpenMP parallelization: true, [anything else]
openmp :=
openmp := true
# MPI parallelization - experimental: true, [anything else]
mpi :=
# Linking: static, shared, dynamic (last two are equal)
......@@ -565,7 +565,7 @@ endif
ifeq ($(lapack),true)
# Mac OS X uses frameworks
ifneq (,$(findstring $(iOS),Darwin))
iLIBS += -framework veclib
iLIBS += -framework Accelerate
else
ifeq ("$(wildcard $(LAPACKDIR)*)","")
$(error Error: LAPACK path '$(LAPACKDIR)' not found.)
......
......@@ -27,11 +27,15 @@ program ED_Kriging
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 ! number of cells
use mo_setVario , only: setVario
nCell, MetSta
use kriging , only: dCS, dS
use mo_setVario , only: setVario, dMatrix
use mo_netcdf , only: NcDataset, NcVariable
use mo_write , only: open_netcdf
use kriging
use mo_EDK , only: EDK
!$ use mo_string_utils, ONLY : num2str
!$ use omp_lib, ONLY : OMP_GET_NUM_THREADS ! OpenMP routines
implicit none
......@@ -42,11 +46,18 @@ program ED_Kriging
integer(i4) :: jDay ! loop variable - current julian day
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
!$OMP PARALLEL
!$ n_threads = OMP_GET_NUM_THREADS()
!$OMP END PARALLEL
!$ print *, 'Run with OpenMP with ', trim(num2str(n_threads)), ' threads.'
call Timer
call ReadDataMain(fname, author_name, vname_data)
......@@ -77,16 +88,27 @@ program ED_Kriging
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 (flagMthTyp)
case (1)
call EDK(jday, iCell)
call EDK(jday, iCell, dCS, MetSta, dS, cell)
case (2)
call OK(jday, iCell)
end select
end do ncellsloop
!$OMP end do
!$OMP end parallel
! correct precipitation values
if (flagVarTyp == 1) then
......
......@@ -13,8 +13,9 @@ subroutine OK(jd,k)
use mainVar
use kriging
use runControl
use varfit, only : beta
use mo_kind, only : i4, dp
use varfit, only : beta
use mo_kind, only : i4, dp
use mo_setVario, only : tVar
! use LSASF_INT
! use lapack95, only: gesv
implicit none
......@@ -24,7 +25,6 @@ subroutine OK(jd,k)
integer(i4) :: ii, jj
integer(i4) :: Nk(nSta), nNmax
real(dp), allocatable :: A (:,:), B(:), X(:)
real(dp) :: tVar
real(dp), allocatable :: lamda(:)
real(dp) :: sumLamda
!
......
......@@ -15,7 +15,7 @@ subroutine OPTI(pmin)
use VarFit
use mo_kind, only: i4, dp
use mo_obj_func, only: obj_func
use mo_nelmin, only: nelmin, nelminrange
use mo_nelmin, only: nelminrange
! parameters for Nelder-Mead algorithm
real(dp) :: pstart(3) ! Starting point for the iteration.
......
module mo_EDK
implicit none
private
public :: EDK
contains
!****************************************************************************
!
! SUBROUTINE: External Drift Kriging
......@@ -11,39 +20,37 @@
! 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)
subroutine EDK(jd, k, dCS, MetSta, dS, cell)
use mo_kind, only : i4, dp
use mainVar
use kriging
use runControl
use mainVar, only : MeteoStation, noDataValue, nSta, thresholdDist
use kriging, only : CellCoarser, dtoS
use runControl, only : flagVarTyp
use varfit, only : beta
! use LSASF_INT
! use mkl95_lapack, only: gesv
! use lapack95, only: gesv
use mo_setVario, only : tVar
implicit none
integer(i4), intent(in) :: jd ! julian day
integer(i4), intent(in) :: k ! cell id
! 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 (:,:), C(:,:), B(:), X(:)
real(dp) :: tVar
real(dp), allocatable :: A (:,:), B(:), X(:) !, C(:,:)
real(dp), allocatable :: lamda(:)
real(dp) :: sumLamda
!
! check DEM
if (nint(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
ll = 0
do i= 1, cell(k)%nNS
do i = 1, cell(k)%nNS
j = cell(k)%listNS(i)
if ( MetSta(j)%z(jd) /= noDataValue ) then
......@@ -55,6 +62,7 @@ subroutine EDK(jd, k)
end if
end do
nNmax = l
!>>>> no value ! avoid indetermination
! avoid 0 value calculations ll == nNmax
! avoid calculations where only 1 station is available
......@@ -136,7 +144,7 @@ subroutine EDK(jd, k)
if (cell(k)%z < 0.0_dp .and. flagVarTyp == 1) then
cell(k)%z = 0.0_dp
end if
deallocate (A,B,X)
deallocate (A,B,X,ipvt)
!
! only one precipiation station available /= 0
else if (ll /= nNmax .AND. nNmax == 1 .AND. flagVarTyp == 1) then
......@@ -158,7 +166,7 @@ subroutine EDK(jd, k)
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
! 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
......@@ -176,171 +184,5 @@ subroutine EDK(jd, k)
end if
!
end subroutine EDK
!
!***********************************************************
!
! tVAR:: Function to calulate variogram at given distance
!
!***********************************************************
real(8) function tVar(h,c0,c,a)
use VarFit, only : vType
use mo_kind, only : dp
real(dp), intent(in) :: h ! distance
real(dp), intent(in) :: c0 ! nugget = beta(1) = XU(1)
real(dp), intent(in) :: c ! sill = beta(2) = XU(2)
real(dp), intent(in) :: a ! range = beta(3) = XU(3)
real(dp) :: r
!
select case (vType)
!
case (1)
! composed: nugget + spherical + sill
r = h/a
if (h == 0.0_dp) then
tVar = 0.0_dp
elseif ( h <= a) then
tVar = c0 + c * (1.5_dp * r - 0.5_dp * r**3)
else
tVar = c0 + c
end if
!
case (2)
! composed: nugget + exponential + sill
r = h/a
tVar = c0 + c * (1.0_dp - dexp(-r))
end select
!
end function tVar
!
!*******************************************************
!
! DMATRIX:: To calculate distance between pairs.......
!
!*******************************************************
subroutine dMatrix
use mo_kind, only : i4, dp
use mainVar
use kriging
use runControl
implicit none
integer(i4) :: i, j, k
integer(i4) :: r, c, ii, jj
integer(i4) :: delta, nTcell
integer(i4) :: NoCellsFiner
real(dp) :: xc, yc
integer(i4), allocatable :: list(:)
!
! Initialize variables
if ( allocated(dCS) ) deallocate (dCS)
if ( allocated(dS) ) deallocate (dS)
if ( allocated(cell)) deallocate (cell)
if ( allocated(dz2S)) deallocate (dz2S)
!
allocate ( dz2S(nSta-1) )
allocate ( dCS(nCell,nSta) )
allocate ( dS(nSta-1) )
allocate ( cell(nCell) )
allocate ( list(nSta) )
!
do i=1,nSta-1
allocate ( dS(i)%S(i+1:nSta) )
allocate ( dz2S(i)%S(i+1:nSta) )
! distance matrix between stations: checked OK
do j=i+1, nSta
dS(i)%S(j) = dsqrt( ( MetSta(i)%x - MetSta(j)%x )**2 + ( MetSta(i)%y - MetSta(j)%y )**2 )
if (dS(i)%S(j) == 0.0_dp) then
print* , 'Stations: ', MetSta(i)%Id, MetSta(j)%Id, ' have the same coordinates, or are repeated. Check LUT.'
!stop
end if
end do
end do
! cell coordinates and elevation : checked OK
! ***************************************
! cell numbering convention (1DIM first)
! c->1 2 ncol
! r
! ---+-----+------...+...-----+
! 1 nr+1
! 2
! ... k
! nr 2nr nCell
! ***************************************
! ii column in finer grid
! yy row in finer grid
! (xc,yc) coordinates of meteogrid
! ***************************************
!
r=1
c=0
NoCellsFiner = (gridMeteo%cellsize / grid%cellsize) **2
xc = gridMeteo%xllcorner + dble(gridMeteo%cellsize) * 0.5_dp
delta = cellFactor / 2
jj = delta
do k=1,nCell
if (r == 1) then
c = c + 1
if (c > 1) then
xc = xc + dble(gridMeteo%cellsize)
jj = jj + cellFactor
end if
yc = gridMeteo%yllcorner + dble(gridMeteo%cellsize) * (dble(gridMeteo%nrows) - 0.5_dp)
ii = delta
else
yc = yc - dble(gridMeteo%cellsize)
ii = ii + cellFactor
end if
cell(k)%x = xc
cell(k)%y = yc
! average of only four DEM cells around centre cell (from lower grid scale upto higher grid cell)
!cell(k)%h = 0.25_dp*(G(ii,jj)%h + G(ii,jj+1)%h + G(ii+1,jj)%h + G(ii+1,jj+1)%h)
!
! average of all DEM cells (from lower grid scale upto higher grid cell)
nTcell = count(G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h > grid%nodata_value )
if (nTcell == 0) then
cell(k)%h = gridMeteo%nodata_value
else
cell(k)%h = sum(G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h, &
G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h /= gridMeteo%nodata_value ) / dble(nTcell)
if ( NoCellsFiner < NoCellsFiner) then
print*, 'Cells with non matching DEM on the finer scale (nodata values in finer grid)'
end if
end if
r=r+1
if (r > gridMeteo%nrows) r = 1
! MZMZMZMZ - delete !!!!!!!
if (cellFactor == 1) then
!print*, ii, jj
cell(k)%h = G(ii+1,jj+1)%h
cycle
end if
! MZMZMZMZ - delete !!!!!!!
end do
! distance matrix cell to stations: checked OK
do j=1, nSta
do i=1,nCell
dCS(i,j) = dsqrt( ( cell(i)%x - MetSta(j)%x )**2 + ( cell(i)%y - MetSta(j)%y )** 2)
end do
end do
! find the closest stations to cell i (any order): checked OK
do i=1,nCell
list = -9
do j=1,nSta
if (dCS(i,j) <= maxDist) list(j) = j
end do
cell(i)%nNS = count(list > -9)
allocate ( cell(i)%listNS( cell(i)%nNS ) )
cell(i)%listNS = pack(list, MASK = list >-9)
end do
!
deallocate (list)
end subroutine dMatrix
end module mo_EDK
......@@ -7,13 +7,13 @@ module mo_obj_func
PUBLIC :: obj_func
contains
function obj_func(p)
use mo_kind, only : i4, dp
use varfit, only : vtype, nbins, gamma, nh
use mo_kind, only : i4, dp
use varfit, only : nbins, gamma, nh
use mo_setVario, only : tVar
implicit none
real(dp), dimension(:), intent(in) :: p
real(dp) :: obj_func
real(dp) :: xu(100)
real(dp) :: gcal, tvar
real(dp) :: gcal
integer(i4) :: k
!
obj_func = 0.0_dp
......
......@@ -90,7 +90,7 @@ module kriging
end type CellCoarser
type(CellCoarser), dimension(:), allocatable :: cell ! EDK output
real(dp), dimension(:,:), allocatable :: dCS ! Euclidean distance between cells -> stations
real(dp), dimension(:,:), allocatable :: dCS ! Euclidean distance between cells -> stations
type dtoS
real(dp), dimension(:), allocatable :: S ! distance to Station j
end type dtoS
......
......@@ -17,6 +17,8 @@ module mo_setVario
private
public :: setVario
public :: dMatrix
public :: tVar
contains
......@@ -57,4 +59,172 @@ subroutine setVario(param)
end if
end subroutine setVario
!
!***********************************************************
!
! tVAR:: Function to calulate variogram at given distance
!
!***********************************************************
real(8) function tVar(h,c0,c,a)
use VarFit, only : vType
use mo_kind, only : dp
real(dp), intent(in) :: h ! distance
real(dp), intent(in) :: c0 ! nugget = beta(1) = XU(1)
real(dp), intent(in) :: c ! sill = beta(2) = XU(2)
real(dp), intent(in) :: a ! range = beta(3) = XU(3)
real(dp) :: r
!
select case (vType)
!
case (1)
! composed: nugget + spherical + sill
r = h/a
if (h == 0.0_dp) then
tVar = 0.0_dp
elseif ( h <= a) then
tVar = c0 + c * (1.5_dp * r - 0.5_dp * r**3)
else
tVar = c0 + c
end if
!
case (2)
! composed: nugget + exponential + sill
r = h/a
tVar = c0 + c * (1.0_dp - dexp(-r))
end select
!
end function tVar
!
!*******************************************************
!
! DMATRIX:: To calculate distance between pairs.......
!
!*******************************************************
subroutine dMatrix
use mo_kind, only : i4, dp
use mainVar
use kriging
use runControl
implicit none
integer(i4) :: i, j, k
integer(i4) :: r, c, ii, jj
integer(i4) :: delta, nTcell
integer(i4) :: NoCellsFiner
real(dp) :: xc, yc
integer(i4), allocatable :: list(:)
!
! Initialize variables
if ( allocated(dCS) ) deallocate (dCS)
if ( allocated(dS) ) deallocate (dS)
if ( allocated(cell)) deallocate (cell)
if ( allocated(dz2S)) deallocate (dz2S)
!
allocate ( dz2S(nSta-1) )
allocate ( dCS(nCell,nSta) )
allocate ( dS(nSta-1) )
allocate ( cell(nCell) )
allocate ( list(nSta) )
!
do i=1,nSta-1
allocate ( dS(i)%S(i+1:nSta) )
allocate ( dz2S(i)%S(i+1:nSta) )
! distance matrix between stations: checked OK
do j=i+1, nSta
dS(i)%S(j) = dsqrt( ( MetSta(i)%x - MetSta(j)%x )**2 + ( MetSta(i)%y - MetSta(j)%y )**2 )
if (dS(i)%S(j) == 0.0_dp) then
print* , 'Stations: ', MetSta(i)%Id, MetSta(j)%Id, ' have the same coordinates, or are repeated. Check LUT.'
!stop
end if
end do
end do
! cell coordinates and elevation : checked OK
! ***************************************
! cell numbering convention (1DIM first)
! c->1 2 ncol
! r
! ---+-----+------...+...-----+
! 1 nr+1
! 2
! ... k
! nr 2nr nCell
! ***************************************
! ii column in finer grid
! yy row in finer grid
! (xc,yc) coordinates of meteogrid
! ***************************************
!
r=1
c=0
NoCellsFiner = (gridMeteo%cellsize / grid%cellsize) **2
xc = gridMeteo%xllcorner + dble(gridMeteo%cellsize) * 0.5_dp
delta = cellFactor / 2
jj = delta
do k=1,nCell
if (r == 1) then
c = c + 1
if (c > 1) then
xc = xc + dble(gridMeteo%cellsize)
jj = jj + cellFactor
end if
yc = gridMeteo%yllcorner + dble(gridMeteo%cellsize) * (dble(gridMeteo%nrows) - 0.5_dp)
ii = delta
else
yc = yc - dble(gridMeteo%cellsize)
ii = ii + cellFactor
end if
cell(k)%x = xc
cell(k)%y = yc
! average of only four DEM cells around centre cell (from lower grid scale upto higher grid cell)
!cell(k)%h = 0.25_dp*(G(ii,jj)%h + G(ii,jj+1)%h + G(ii+1,jj)%h + G(ii+1,jj+1)%h)
!
! average of all DEM cells (from lower grid scale upto higher grid cell)
nTcell = count(G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h > grid%nodata_value )
if (nTcell == 0) then
cell(k)%h = gridMeteo%nodata_value
else
cell(k)%h = sum(G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h, &
G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h /= gridMeteo%nodata_value ) / dble(nTcell)
if ( NoCellsFiner < NoCellsFiner) then
print*, 'Cells with non matching DEM on the finer scale (nodata values in finer grid)'
end if
end if
r=r+1
if (r > gridMeteo%nrows) r = 1
! MZMZMZMZ - delete !!!!!!!
if (cellFactor == 1) then
!print*, ii, jj