Commit 99be6730 authored by Friedrich Boeing's avatar Friedrich Boeing

upload old edk version using imsl server

parent ee30d840
File added
# This cmake file is not meant to be edited for a special setup.
# For special setups use cache line files or command line options, as described a few
# lines ahead concerning module independend builds
cmake_minimum_required(VERSION 3.5)
project(edk Fortran)
# 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:
# $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")
message(STATUS "build independend of module system ${CMAKE_BUILD_MODULE_SYSTEM_INDEPENDEND}")
# set specific place where to search for the netCDF directory
set(CMAKE_NETCDF_DIR " " CACHE STRING "set set specific place where to search for the netCDF directory")
message(STATUS "search in additional directory ${CMAKE_NETCDF_DIR} for netCDF")
# The variable "CMAKE_WITH_MPI" can be set before executing cmake via a cache command:
# $cmake -DCMAKE_WITH_MPI:STRING=ON ..
# or in a cache file:
# $cmake -C ../CMakeCacheFiles/example
# or after executing CMake editing the CMakeCache.txt, preferably with a corresponding cmake editor i.e. ccmake
set(CMAKE_WITH_MPI OFF CACHE STRING "build the module with MPI, so it can be executed using mpirun")
# same with OpenMP
set(CMAKE_WITH_OpenMP OFF CACHE STRING "build the module with OpenMP parallelization")
# same with lapack
set(CMAKE_WITH_LAPACK OFF CACHE STRING "build the module with lapack library")
# additional cmake-modules created for the purpose of finding netCDF or other libraries ly in the source_directory in
# a folder named cmake-modules. This command tells cmake to search there for Find<module>.cmake files
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake-modules)
set (NETCDF_F90 "YES")
# the FindNetCDFF.cmake file can be found after we added the cmake-modules folder to the CMAKE_MODULE_PATH
# the build fails, if it is not present
find_package(NetCDFF REQUIRED)
# from that module we gain the following variables:
# NETCDF_INCLUDES : the include directory
# NETCDF_LINK_LIBRARIES : the absolute path to and with the libraries
# NETCDF_CFLAGS_OTHER : additional compilation flags
# NETCDF_LDFLAGS_OTHER : additional linking flags
# if cmake provides a findLIBRARY module, this gets invoked via find_package(LIBRARY)
if (CMAKE_WITH_MPI)
# find if there is an MPI setup on the system and if so, set corresponding variables
find_package(MPI)
if (NOT ${MPI_Fortran_FOUND})
message(FATAL_ERROR "MPI required but not found")
else()
message(STATUS "found MPI_Fortran_COMPILER ${MPI_Fortran_COMPILER}")
endif()
endif()
if (CMAKE_WITH_OpenMP)
# find if there is an OpenMP setup on the system and if so, set corresponding variables
find_package(OpenMP)
if (NOT ${OpenMP_Fortran_FOUND})
message(FATAL_ERROR "OpenMP required but not found")
endif()
endif()
if (CMAKE_WITH_LAPACK)
# find if there is an LAPACK library on the system and if so, set corresponding variables
find_package(LAPACK)
if (NOT ${LAPACK_FOUND})
message(FATAL_ERROR "lapack required but not found")
endif()
endif()
include_directories(${NETCDF_INCLUDES} ${MPI_Fortran_INCLUDE_PATH} ${OpenMP_Fortran_LIBRARY})
# ifort and gfortran need the flag -cpp to interpret definitions like -DMRM2MHM
# the nag compiler is not able to interpret the flag -cpp but can interpret these definitions anyway
# so we check whether the compiler is able to use the flag -cpp
# for that we need the module CheckFortranCompilerFlag
include(CheckFortranCompilerFlag)
CHECK_Fortran_COMPILER_FLAG("-cpp" CPP_FLAG)
# if the flag exists, we add it to the compilation flags
if (CPP_FLAG)
set(ADDITIONAL_GCC_FLAGS "-cpp")
endif()
# this function adds definitions but also creates a corresponding CMAKE variable with CACHE STRING
# i.e.:
# The variable "${defCMakeName}" can be set before executing cmake via a cache command cmake -D...
# or in a cache file:
# $cmake -C ../CMakeCacheFiles/example
# or after executing CMake editing the CMakeCache.txt, preferably with a corresponding cmake editor i.e. ccmake
# cmake ..
function(cpp_definitions defName defCMakeName value cacheString)
set(${defCMakeName} "${value}" CACHE STRING "${cacheString}")
if (${defCMakeName})
add_definitions("${defName}")
endif()
endfunction()
# Add definitions. These should later be set via the cache line file and only
# have a default value here. These part only concerns user build definitions like MRM2MHM
#add_definitions(-DMRM2MHM=.true.)
cpp_definitions("-DABSOFT" "CMAKE_ABSOFT" "OFF" "Documentation to be added. If you you are developer, you might edit this string in CMakeLists.txt")
# Compile.
#
# This is not recommended but wished by some members of our working group. With this command
# all files in src with the ending f90 or h are written into sources, and therefore linked
# together to the executable, if relevant or not
file(GLOB_RECURSE sources src/*.f90 src/*.h)
# this command is able to create dependencies, compile and add the sources in the right order
add_executable(edk ${sources})
# the libraries are added to the executable by the linker, using the full path and using the
# rpath option, except the libraries are located in ${CMAKE_Fortran_IMPLICIT_LINK_DIRECTORIES}.
target_link_libraries(edk ${NETCDF_LINK_LIBRARIES} ${MPI_Fortran_LIBRARIES} ${OpenMP_Fortran_LIBRARIES} ${LAPACK_LIBRARIES})
set_property(TARGET edk PROPERTY COMPILE_FLAGS "${CMAKE_Fortran_FLAGS} ${ADDITIONAL_GCC_FLAGS} ${NETCDF_CFLAGS_OTHER} ${MPI_Fortran_COMPILE_FLAGS} ${OpenMP_Fortran_FLAGS}")
set_property(TARGET edk PROPERTY LINK_FLAGS "${NETCDF_LDFLAGS_OTHER} ${MPI_Fortran_LINK_FLAGS} ${OpenMP_Fortran_FLAGS} ${LAPACK_LINKER_FLAGS}")
# set compiling flags for debug and realese version
if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-form -ffixed-line-length-132")
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -pedantic-errors -Wall -W -O -g -Wno-maybe-uninitialized")
set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -O3")
cpp_definitions("-DGFORTRAN" "CMAKE_GFORTRAN" "ON" "Code exchange for gfortran compiler dependent issues")
endif()
if(CMAKE_Fortran_COMPILER_ID MATCHES "Intel")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -nofixed -assume byterecl -fp-model source -m64 -assume realloc-lhs ") # precise -> source: suppress warning, computation identical
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -warn all -g -debug -traceback -fp-stack-check -O0 -debug -check all")
set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -O3 -qoverride-limits")
endif()
message(STATUS "the following debug flags will be used: ${CMAKE_Fortran_FLAGS_DEBUG}")
# Usually that works fine, except, one is on a module system and tries to execute the executable
# in the end without having the modules loaded. A workaround is provided using the variable
# CMAKE_BUILD_MODULE_SYSTEM_INDEPENDEND
# if this variable is set to ON (do not set the variable inside of this cmake file), then the
# paths are added to the INSTALL_RPATH, and via the second command also to the build.
# It is a bit of a mess and workaround though.
if (CMAKE_BUILD_MODULE_SYSTEM_INDEPENDEND)
set_target_properties(edk PROPERTIES INSTALL_RPATH "${CMAKE_Fortran_IMPLICIT_LINK_DIRECTORIES}")
set_target_properties(edk PROPERTIES BUILD_WITH_INSTALL_RPATH ON)
endif()
#set_target_properties(edk PROPERTIES SKIP_BUILD_RPATH OFF)
#set_target_properties(edk PROPERTIES INSTALL_RPATH_USE_LINKPATH ON)
# possibility to create symlinks to the build or to the source directory, so a copy of mhm does not have to be done. On the other
# hand, only for the test basins the executable would not be copied or linked to a predictable location. So it may be a good
# idea to not copy it with the cmake setup
#add_custom_target(link_namelists ALL "${CMAKE_COMMAND}" -E create_symlink "${CMAKE_CURRENT_SOURCE_DIR}/mhm.nml" "${CMAKE_CURRENT_BINARY_DIR}/mhm.nml")
......@@ -4,12 +4,12 @@
! http://www.mohid.com/wiki/index.php?title=PROJ4
subroutine CoorTrans(lon, lat, x, y)
! use proj4
use proj4
implicit none
real(8), intent(in) :: lon, lat
real(8), intent(out) :: x, y
integer(4) :: status
! type(prj90_projection) :: proj
type(prj90_projection) :: proj
character(len=20), dimension(8) :: params
!
! define parameters of the aiming coordinate system
......@@ -23,32 +23,28 @@ subroutine CoorTrans(lon, lat, x, y)
params(7) = 'x_0=4500000' ! false easting
params(8) = 'y_0=0' ! false northing
!
! dummy variables
x = 0.
y = 0.
status=prj90_init(proj,params)
if (status.ne.PRJ90_NOERR) then
write(*,*) prj90_strerrno(status)
stop
end if
!
! status=prj90_init(proj,params)
! if (status.ne.PRJ90_NOERR) then
! write(*,*) prj90_strerrno(status)
! stop
! end if
! !
! ! lat,lon --> GK4
! status = prj90_fwd(proj, lon, lat, x, y)
! if (status.ne.PRJ90_NOERR) then
! write(*,*) prj90_strerrno(status)
! stop
! end if
! lat,lon --> GK4
status = prj90_fwd(proj, lon, lat, x, y)
if (status.ne.PRJ90_NOERR) then
write(*,*) prj90_strerrno(status)
stop
end if
end subroutine CoorTrans
subroutine CoorTransInv(x, y, lon, lat)
! use proj4
use proj4
implicit none
real(8), intent(in) :: x, y
real(8), intent(out) :: lon, lat
integer(4) :: status
! type(prj90_projection) :: proj
type(prj90_projection) :: proj
character(len=20), dimension(8) :: params
!
! define parameters of the aiming coordinate system
......@@ -62,26 +58,22 @@ subroutine CoorTransInv(x, y, lon, lat)
params(7) = 'x_0=4500000' ! false easting
params(8) = 'y_0=0' ! false northing
!
! dummy variables
lon = 0.
lat = 0.
!print*, ' '
!print*, 'CoordTransInv: before init'
status=prj90_init(proj,params)
!print*, 'CoordTransInv: after init'
if (status.ne.PRJ90_NOERR) then
write(*,*) prj90_strerrno(status)
stop
end if
!
! !print*, ' '
! !print*, 'CoordTransInv: before init'
! status=prj90_init(proj,params)
! !print*, 'CoordTransInv: after init'
! if (status.ne.PRJ90_NOERR) then
! write(*,*) prj90_strerrno(status)
! stop
! end if
! !
! ! GK 4 --> lat, lon
! !print*, 'CoordTransInv: before transformation', x, y
! status = prj90_inv(proj, x, y, lon, lat)
! !print*, 'CoordTransInv: after transformation', lon, lat
! !print*, ' '
! if (status.ne.PRJ90_NOERR) then
! write(*,*) prj90_strerrno(status)
! stop
! end if
! GK 4 --> lat, lon
!print*, 'CoordTransInv: before transformation', x, y
status = prj90_inv(proj, x, y, lon, lat)
!print*, 'CoordTransInv: after transformation', lon, lat
!print*, ' '
if (status.ne.PRJ90_NOERR) then
write(*,*) prj90_strerrno(status)
stop
end if
end subroutine CoorTransInv
!**********************************************************************************
! VARIOGRAM: Seting or estimating and fitting
! PURPOSE:
! 1) Set variagram from DB for a block, or
! 2.1) Estimate an empirical semi variogram for daily met. values
! 2.2) Fitting a teoretical variogram
! 2.3) Set variogram for a block
! WHERE:
! UPDATES:
! Created Sa 19.02.2004 main structure
! Last update 12.04.2006
!**********************************************************************************
module mo_setVario
!****************************************************************************
!
! 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)
use mo_kind, only : i4, dp
use mainVar
use kriging
use runControl
use varfit, only : beta
use LSASF_INT
implicit none
integer(i4), intent(in) :: jd ! julian 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) :: tVar
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
ll = 0
private
do i= 1, cell(k)%nNS
j = cell(k)%listNS(i)
if ( MetSta(j)%z(jd) /= noDataValue ) then
public :: setVario
public :: dMatrix
public :: tVar
!*** 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
!
! initialize matrices
allocate (A(nNmax+2,nNmax+2), B(nNmax+2), X(nNmax+2))
!
! 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
A(i,nNmax+2) = MetSta(ii)%h
B(i)=tVar(dCS(k,ii),beta(1),beta(2),beta(3))
end do
A(nNmax,nNmax+1) = 1.0_dp
ii=Nk(nNmax)
A(nNmax,nNmax+2) = MetSta(ii)%h
B(nNmax)=tVar(dCS(k,ii),beta(1),beta(2),beta(3))
!
B(nNmax+1) = 1.0_dp
B(nNmax+2) = cell(k)%h
contains
!NOTE: only the upper triangular matrix is needed!
call D_LSASF (A, B, X)
!
! The BLUE of z is then:
cell(k)%z = 0.
do i=1,nNmax
ii=Nk(i)
cell(k)%z = cell(k)%z + 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 (cell(k)%z < 0.0_dp .and. flagVarTyp == 1) then
cell(k)%z = 0.0_dp
end if
deallocate (A,B,X)
!
! only one precipiation station available /= 0
else if (ll /= nNmax .AND. nNmax == 1 .AND. flagVarTyp == 1) then
ii = Nk(1)
! problably convective rain at the only station
if (dCS(k,ii) <= thresholdDist) then
cell(k)%z = MetSta(ii)%z(jd)
else
cell(k)%z = 0.0_dp
end if
subroutine setVario(param)
use runControl
use VarFit
use mainVar
use mo_EmpVar, only: EmpVar
use mo_kind, only: i4, dp
implicit none
real(dp), intent(out) :: param(3)
integer(i4) :: jd, y, y0
!
! Estimation
if (flagVario) then
! Initialize
nobs=0 ; m0=0.0_dp ; v0=0.0_dp
y0 = 0
!
do jd= jStart, jEnd
y = floor(float(jd-jStart)/365.)
if (y > y0 ) then
y0 = y
print*, 'VarFit. Processing ', yStart+y0-1 , '...'
end if
! update Variogram bins daily
call EmpVar(jd, .true.)
end do
!
! variance
v0 = v0 / real(nobs, dp)
! optimize parameters
! open(UNIT = 6,FILE = 'Report_OPT.sol', STATUS='REPLACE')
! only one station available, which has values /= 0
else if (ll /= nNmax .AND. nNmax == 1 .AND. flagVarTyp /= 1) then
ii = Nk(1)
cell(k)%z = MetSta(ii)%z(jd)
print *, "ST: replace old GRG2 opti with something better"
call opti(param)
! 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 setVario
!
end subroutine EDK
!
!***********************************************************
!
......@@ -80,7 +161,7 @@ real(8) function tVar(h,c0,c,a)
! composed: nugget + spherical + sill
r = h/a
if (h == 0.0_dp) then
tVar = c0 ! 0.0_dp