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

Checking Lapack version

- added dgesv for Mac OS (untested)
- substituted some dfloat with real( , dp) commands
parent 724358f9
......@@ -19,7 +19,7 @@ subroutine EDK(jd,k)
use varfit, only : beta
! use LSASF_INT
! use mkl95_lapack, only: gesv
use lapack95, only: gesv
! use lapack95, only: gesv
implicit none
integer(i4), intent(in) :: jd ! julian day
......@@ -104,7 +104,8 @@ subroutine EDK(jd,k)
! NOTE: only the upper triangular matrix is needed!
! call D_LSASF (A, B, X)
call gesv(A, B)
! call gesv(A, B) ! EVE CentOS 7
call dgesv(A, B) ! MacOS Sierra: Accelerate Framework
X = B
!
! The BLUE of z is then:
......
......@@ -58,7 +58,7 @@ program ED_Kriging
! estimate variogram
call setVario
! write variogram
if (flagVario == .TRUE.) call WriteDataMeteo(0,0,2)
if (flagVario) call WriteDataMeteo(0,0,2)
!
if (flagEDK) then
! open netcdf if necessary
......
......@@ -46,7 +46,7 @@ subroutine EmpVar(jd, flagMax)
if ( MetSta(i)%z(jd) == noDataValue ) cycle
nobs = nobs + 1
delta = MetSta(i)%z(jd) - m0
m0 = m0 + delta / dfloat(nobs)
m0 = m0 + delta / real(nobs, dp)
v0 = v0 + delta * (MetSta(i)%z(jd) - m0)
!write (999,'(4i6, 3f15.4)') y,d,i, nobs, MetSta(i)%z(jd), m0, v0
end do
......@@ -125,12 +125,12 @@ subroutine EmpVar(jd, flagMax)
if (Nh(k) > 0) then
ni=ni+1
! Classsical
gamma(k,2)=gamma(k,2)/2.0_dp/dfloat(Nh(k))
gamma(k,2)=gamma(k,2)/2.0_dp/real(Nh(k), dp)
!
! 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)=dfloat(k)*dh-dfloat(ni)*dh/2.0_dp
gamma(k,1)=real(k, dp)*dh-real(ni, dp)*dh/2.0_dp
ni=0
else
ni=ni+1
......
......@@ -16,7 +16,7 @@ subroutine OK(jd,k)
use varfit, only : beta
use mo_kind, only : i4, dp
! use LSASF_INT
use lapack95, only: gesv
! use lapack95, only: gesv
implicit none
integer(i4), intent(in) :: jd ! day
integer(i4), intent(in) :: k ! cell id
......
......@@ -84,9 +84,9 @@ module kriging
integer(i4) :: nNS ! No. Nearest Stations (NS) d <= maxDist
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) :: z ! z values to be interpolated (OUTPUT)
real(dp) :: y ! y- coordinate
real(dp) :: h ! (estimated) elevation [m] (from the nearest cells DEM)
real(sp) :: z ! z values to be interpolated (OUTPUT)
end type CellCoarser
type(CellCoarser), dimension(:), allocatable :: cell ! EDK output
......
!**********************************************************************************
! 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
!**********************************************************************************
subroutine setVario
use runControl
use VarFit
use mainVar
use mo_kind, only: i4, dp
implicit none
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 / dfloat(nobs)
! optimize parameters
open(UNIT = 6,FILE = 'Report_OPT.sol', STATUS='REPLACE')
call opti
end if
end subroutine setVario
!**********************************************************************************
! 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
!**********************************************************************************
subroutine setVario
use runControl
use VarFit
use mainVar
use mo_kind, only: i4, dp
implicit none
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')
print *, "ST: replace old GRG2 opti with something better"
! call opti
! fit a theoretical variogram
end if
end subroutine setVario
......@@ -29,12 +29,12 @@ subroutine stats
end do
!
error = zObs-zCal
if ( ne > 1 ) then
zObsMean = sum(zObs)/dfloat(ne)
zCalMean = sum(zCal)/dfloat(ne)
if ( ne > 1 ) then
zObsMean = sum(zObs)/real(ne, dp)
zCalMean = sum(zCal)/real(ne, dp)
sumP = dot_product(zObs,zCal)
zObsVar = dot_product(zObs,zObs) - dfloat(ne) * zObsMean * zObsMean
zCalVar = dot_product(zCal,zCal) - dfloat(ne) * zCalMean * zCalMean
zObsVar = dot_product(zObs,zObs) - real(ne, dp) * zObsMean * zObsMean
zCalVar = dot_product(zCal,zCal) - real(ne, dp) * zCalMean * zCalMean
SSE = dot_product(error,error)
denom = zObs - zObsMean
NSE_denom = dot_product(denom,denom)
......@@ -53,7 +53,7 @@ subroutine stats
E(1) = zCalMean - zObsMean
!
! MSE
E(2) = SSE/dfloat(ne)
E(2) = SSE/real(ne, dp)
!
! RMSE
if ( E(2) > 0.0_dp ) then
......@@ -71,16 +71,16 @@ subroutine stats
!
! MAE
E(5)= sum(abs(error))
E(5)= E(5)/DFLOAT(ne)
E(5)= E(5)/real(ne, dp)
!
! RMAE
E(6)=E(5)/zObsMean
!
! r
E(7)= (sumP-dfloat(ne) * zCalMean * zObsMean) / dsqrt(zCalVar * zObsVar)
!
! NSE
E(8)= 1.0_dp - (SSE/NSE_denom)
E(7)= (sumP-real(ne, dp) * zCalMean * zObsMean) / dsqrt(zCalVar * zObsVar)
!
! NSE
E(8)= 1.0_dp - (SSE/NSE_denom)
else
E = small
end if
......
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