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

first implementation of openMP parallelization

parent 61f99358
......@@ -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.)
......
......@@ -29,11 +29,13 @@ program ED_Kriging
grid, gridMeteo, & ! grid properties of input and output grid
nCell, MetSta
use kriging , only: dCS, dS
use mo_setVario , only: setVario
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, dMatrix
use mo_EDK , only: EDK
!$ use mo_string_utils, ONLY : num2str
!$ use omp_lib, ONLY : OMP_GET_NUM_THREADS ! OpenMP routines
implicit none
......@@ -44,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)
......@@ -79,8 +88,11 @@ 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
......@@ -95,6 +107,8 @@ program ED_Kriging
call OK(jday, iCell)
end select
end do ncellsloop
!$OMP end do
!$OMP end parallel
! correct precipitation values
if (flagVarTyp == 1) then
......
......@@ -13,9 +13,9 @@ subroutine OK(jd,k)
use mainVar
use kriging
use runControl
use varfit, only : beta
use mo_kind, only : i4, dp
use mo_EDK, only : tVar
use varfit, only : beta
use mo_kind, only : i4, dp
use mo_setVario, only : tVar
! use LSASF_INT
! use lapack95, only: gesv
implicit none
......
......@@ -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.
......
......@@ -5,8 +5,6 @@ module mo_EDK
private
public :: EDK
public :: dMatrix
public :: tVar
contains
!****************************************************************************
......@@ -28,13 +26,13 @@ subroutine EDK(jd, k, dCS, MetSta, dS, cell)
use kriging, only : CellCoarser, dtoS
use runControl, only : flagVarTyp
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), dimension(:, :), intent(in) :: dCS ! distances matrix
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
......@@ -44,7 +42,7 @@ subroutine EDK(jd, k, dCS, MetSta, dS, cell)
integer(i4) :: ii, jj
integer(i4) :: Nk(nSta), nNmax
integer(i4), allocatable :: ipvt(:)
real(dp), allocatable :: A (:,:), C(:,:), B(:), X(:)
real(dp), allocatable :: A (:,:), B(:), X(:) !, C(:,:)
real(dp), allocatable :: lamda(:)
real(dp) :: sumLamda
!
......@@ -186,172 +184,5 @@ subroutine EDK(jd, k, dCS, MetSta, dS, cell)
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,12 @@ 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_EDK, only : tVar
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
integer(i4) :: k
!
......
......@@ -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
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_setVario
......@@ -4,9 +4,9 @@
!
! **************************************************************************
subroutine stats
use varFit, only : E, beta, gamma, nbins
use varFit, only : E, beta, gamma
use mo_kind, only : i4, dp
use mo_EDK, only : tVar
use mo_setVario, only : tVar
implicit none
integer(i4), parameter :: incx = 1
integer(i4) :: k
......
......@@ -18,12 +18,12 @@
! Last Update Sa
!**************************************************************************
subroutine WriteDataMeteo(y,d,wFlag)
use mo_kind, only : i4, sp, dp
use mo_kind, only : i4, dp
use mainVar
use runControl
use kriging
use VarFit
use mo_EDK, only : tvar
use mo_setVario, only : tvar
!
implicit none
integer(i4), intent (in) :: y, d, wFlag
......
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