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

Commit 61f99358 authored by Stephan Thober's avatar Stephan Thober
Browse files

modularized EDK subroutine

parent 6a932321
......@@ -27,11 +27,13 @@ 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
nCell, MetSta
use kriging , only: dCS, dS
use mo_setVario , only: setVario
use mo_netcdf , only: NcDataset, NcVariable
use mo_write , only: open_netcdf
use kriging
use mo_EDK , only: EDK, dMatrix
implicit none
......@@ -78,10 +80,16 @@ program ED_Kriging
print *, 'YEAR: ',year, 'DOY: ', doy
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)
......
......@@ -15,6 +15,7 @@ subroutine OK(jd,k)
use runControl
use varfit, only : beta
use mo_kind, only : i4, dp
use mo_EDK, 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
!
......
module mo_EDK
implicit none
private
public :: EDK
public :: dMatrix
public :: tVar
contains
!****************************************************************************
!
! SUBROUTINE: External Drift Kriging
......@@ -11,39 +22,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
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), dimension(:, :), 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 :: 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 +64,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 +146,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 +168,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
......@@ -344,3 +354,4 @@ subroutine dMatrix
end subroutine dMatrix
end module mo_EDK
......@@ -9,11 +9,12 @@ contains
function obj_func(p)
use mo_kind, only : i4, dp
use varfit, only : vtype, nbins, gamma, nh
use mo_EDK, 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
......
......@@ -6,6 +6,7 @@
subroutine stats
use varFit, only : E, beta, gamma, nbins
use mo_kind, only : i4, dp
use mo_EDK, only : tVar
implicit none
integer(i4), parameter :: incx = 1
integer(i4) :: k
......@@ -16,7 +17,6 @@ subroutine stats
real(dp), dimension(:), allocatable :: error, denom
real(dp), dimension(:), allocatable :: zCal, zObs
real(dp), parameter :: small = -9.999d3
real(dp) :: tVar
!
!Initialize
......
......@@ -23,6 +23,7 @@ subroutine WriteDataMeteo(y,d,wFlag)
use runControl
use kriging
use VarFit
use mo_EDK, only : tvar
!
implicit none
integer(i4), intent (in) :: y, d, wFlag
......@@ -31,7 +32,6 @@ subroutine WriteDataMeteo(y,d,wFlag)
character(256) :: dummy
character(256) :: fileName
logical :: wasOpened
real(dp) :: tVar
!
select case (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