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

Commit 2ecdb269 authored by Stephan Thober's avatar Stephan Thober
Browse files

modularized Reading

- changed name of namelist file from main.dat to edk.nml
- updated namelist:
	- included flag for correcting negative interpolation values
	- included flag for setting distant interpolation values to zero
	- included variable for interpolation method, 0 -> no interpolation
	- removed flagVarTyp and flagEDK
parent 0cd30f01
......@@ -21,9 +21,8 @@ program ED_Kriging
use mo_kind , only: i4, dp
use mo_julian , only: NDAYS, NDYIN
use runControl , only: outputformat, & ! outputformat either 'nc' or 'bin'
flagEDK, flagMthTyp, & ! flag for activate kriging, flag for 'OK' or 'EDK'
flagVarTyp, & ! pre or temp
use runControl , only: flagEDK, interMth, & ! flag for activate kriging, flag for 'OK' or 'EDK'
correctNeg, & ! pre or temp
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
......@@ -71,7 +70,7 @@ program ED_Kriging
call message('')
call message(' >>> Reading data')
call message('')
call ReadDataMain(fname, author_name, vname_data)
call ReadDataMain
! read DEM
call ReadDEM
......@@ -113,13 +112,13 @@ program ED_Kriging
if (flagEDK) then
if (interMth .gt. 0) then
itimer = 4
call timer_start(itimer)
call message(' >>> Perform interpolation')
call message('')
! open netcdf if necessary
call open_netcdf(fname, author_name, vname_data, gridMeteo%ncols, gridMeteo%nrows, yStart, nc_out, nc_data, nc_time)
call open_netcdf(nc_out, nc_data, nc_time)
timeloop: do jday = jStart, jEnd
......@@ -139,19 +138,18 @@ program ED_Kriging
cycle
end if
! interploation
select case (flagMthTyp)
select case (interMth)
case (1)
call EDK(jday, iCell, dCS, MetSta, dS, cell)
case (2)
call OK(jday, iCell)
case (2)
call EDK(jday, iCell, dCS, MetSta, dS, cell)
end select
end do ncellsloop
!$OMP end do
!$OMP end parallel
! correct precipitation values
if (flagVarTyp == 1) then
! correct negative values
if (correctNeg) then
where ((cell(:)%z .LT. 0.0_sp) .AND. (cell(:)%z .GT. real(grid%nodata_value, sp)) )
cell(:)%z = 0.0_sp
end where
......
......@@ -26,7 +26,6 @@ subroutine EmpVar(jd, flagMax)
use mo_kind , only : i4, dp
use kriging
use VarFit
use runControl, only : flagVarTyp
implicit none
integer(i4), intent(in) :: jd
logical, intent(in) :: flagMax
......@@ -71,18 +70,19 @@ subroutine EmpVar(jd, flagMax)
gamma = 0.0_dp
end if
!
print *, '***WARNING: Removal of outliers in the estimation of the variogram is deactivated'
do i=1,nSta-1
do j=i+1,nSta
if (dz2S(i)%S(j) /= noDataValue ) then
! take values up to max distance
if (dS(i)%S(j) > hMax ) cycle
! remove outliers for the estimation of the variogram
if (flagVarTyp == 2 ) then
if ( dabs( MetSta(i)%h - MetSta(j)%h ) / dS(i)%S(j) > gradE ) then
! write(999,*), 'pair removed', MetSta(i)%id, MetSta(j)%id
cycle
end if
end if
! ! remove outliers for the estimation of the variogram
! if (flagVarTyp == 2 ) then
! if ( dabs( MetSta(i)%h - MetSta(j)%h ) / dS(i)%S(j) > gradE ) then
! ! write(999,*), 'pair removed', MetSta(i)%id, MetSta(j)%id
! cycle
! end if
! end if
!
k=max(1, ceiling(dS(i)%S(j)/dh))
Nh(k)=Nh(k)+1
......
......@@ -85,21 +85,18 @@ subroutine OK(jd,k)
ii=Nk(i)
cell(k)%z = cell(k)%z + real(X(i) * MetSta(ii)%z(jd))
end do
! 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
! 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)
!
! only one station available /= 0
else if (ll /= nNmax .AND. nNmax == 1) then
ii = Nk(1)
cell(k)%z = MetSta(ii)%z(jd)
!
! for precipitation, distant values are set to zero
if (distZero .and. dCS(k,ii) .gt. thresholdDist) then
cell(k)%z = 0.0_dp
end if
!
! if all stations have the value 0
! if all stations have the value 0
else if (ll == nNmax ) then
......
......@@ -7,21 +7,21 @@
! Created Sa 21.03.2006
! Last Update Sa 11.06.2010
!**************************************************************************
subroutine ReadDataMain(FileOut, author_name, variable_name)
use mainVar
subroutine ReadDataMain()
use mainVar, only: noDataValue, cellFactor, DataConvertFactor, &
yStart, mStart, dStart, yEnd, mEnd, dEnd
use mo_kind, only: i4, dp
use mo_julian , only: NDAYS
use kriging
use VarFit
use runControl
use NetCDFVar, only: variable_unit, variable_long_name
use kriging, only: maxDist
use VarFit, only: vType, nParam, dh, hMax, beta
use runControl, only: interMth, fnameDEM, DataPathOut, DataPathIn, fNameSta, correctNeg, &
distZero, flagVario, fNameVario, flagEDK
use NetCDFVar, only: FileOut, author_name, variable_name, variable_unit, variable_long_name
use mo_message, only: message
implicit none
!
character(256), intent(out) :: FileOut
character(256), intent(out) :: author_name
character(256), intent(out) :: variable_name
integer(i4) :: i, ios
character(256) :: dummy, fileName
!
......@@ -29,18 +29,33 @@ subroutine ReadDataMain(FileOut, author_name, variable_name)
! namelist definition
!===============================================================
!
namelist/mainVars/flagMthTyp, flagVarTyp, noDataValue, DataPathIn, fNameDEM, &
DataPathOut, FileOut, fNameSTA, cellFactor, outputformat, DataConvertFactor, flagEDK, &
author_name, variable_name, variable_unit, variable_long_name, &
namelist/mainVars/noDataValue, DataPathIn, fNameDEM, &
DataPathOut, FileOut, fNameSTA, cellFactor, DataConvertFactor, InterMth, correctNeg, &
distZero, author_name, variable_name, variable_unit, variable_long_name, &
yStart, mStart, dStart, yEnd, mEnd, dEnd, maxDist, flagVario, vType, nParam, &
fNameVario, dh, hMax, author_name, variable_name
fNameVario, dh, hMax
!
! -----------------------------------------------------------------------
! MAIN.DAT
! -----------------------------------------------------------------------
open(unit=10, file='main.dat', STATUS='OLD', ACTION='read')
open(unit=10, file='edk.nml', STATUS='OLD', ACTION='read')
read(10, nml=mainVars)
close(10)
! -----------------------------------------------------------------------
! consistency checks
! -----------------------------------------------------------------------
if ((interMth .gt. 2) .or. (interMth .lt. 0)) then
call message(' >>> Interpolated Method Value not valid')
call message(' Please choose one of the following:')
call message(' 0 - no interpolation')
call message(' 1 - ordinary kriging')
call message(' 2 - external drift kriging')
call message('')
stop 1
end if
FileOut = trim(DataPathOut) // trim(FileOut)
!
! *************************************
......
......@@ -24,7 +24,7 @@ subroutine EDK(jd, k, dCS, MetSta, dS, cell)
use mo_kind, only : i4, dp
use mainVar, only : MeteoStation, noDataValue, nSta, thresholdDist
use kriging, only : CellCoarser, dtoS
use runControl, only : flagVarTyp
use runControl, only : correctNeg, distZero
use varfit, only : beta
use mo_setVario, only : tVar
implicit none
......@@ -141,27 +141,22 @@ subroutine EDK(jd, k, dCS, MetSta, dS, cell)
! 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
if (correctNeg) then
cell(k)%z = 0.0_dp
end if
deallocate (A,B,X,ipvt)
!
! 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
! 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)
! if all stations have the value 0
! only one station available /= 0
else if (ll /= nNmax .AND. nNmax == 1) then
ii = Nk(1)
cell(k)%z = MetSta(ii)%z(jd)
!
! for precipitation, distant values are set to zero
if (distZero .and. dCS(k,ii) .gt. thresholdDist) then
cell(k)%z = 0.0_dp
end if
!
! if all stations have the value 0
else if (ll == nNmax ) then
cell(k)%z = 0.0_dp
......
......@@ -8,20 +8,16 @@ module mo_write
CONTAINS
subroutine open_netcdf(fname, author_name, vname_data, ncols, nrows, year_start, nc, var_data, var_time)
subroutine open_netcdf(nc, var_data, var_time)
use mo_kind, only: i4, sp, dp
use mo_netcdf, only: NcDataset, NcDimension, NcVariable
use mo_string_utils, only: num2str
use mainVar, only: gridMeteo, yStart
use NetCDFVar, only: fileOut, author_name, variable_name
implicit none
character(256), intent(in) :: fname
character(256), intent(in) :: author_name
character(256), intent(in) :: vname_data
integer(i4), intent(in) :: ncols
integer(i4), intent(in) :: nrows
integer(i4), intent(in) :: year_start
type(NcDataset), intent(out) :: nc
type(NcVariable), intent(out) :: var_time, var_data
......@@ -29,32 +25,26 @@ CONTAINS
type(NcVariable) :: var_lon, var_lat
! 1.1 create a file
nc = NcDataset(trim(fname), "w")
nc = NcDataset(trim(fileOut), "w")
! create dimensions
dim_x = nc%setDimension("x", ncols)
dim_y = nc%setDimension("y", nrows)
dim_x = nc%setDimension("x", gridMeteo%ncols)
dim_y = nc%setDimension("y", gridMeteo%nrows)
dim_time = nc%setDimension("time", -1)
! create variables
var_time = nc%setVariable('time', "i32", (/dim_time/))
! var_lat = nc%setVariable(vname_lat, "f32", (/dim_x, dim_y/))
! var_lon = nc%setVariable(vname_lon , "f32", (/dim_x, dim_y/))
var_data = nc%setVariable(vname_data, "f64", (/dim_y, dim_x, dim_time/))
var_data = nc%setVariable(variable_name, "f64", (/dim_y, dim_x, dim_time/))
! add some variable attributes
call var_time%setAttribute("units", "days since " // trim(num2str(year_start - 1, form='(I4)')) // "-12-31 12:00:00")
call var_time%setAttribute("units", "days since " // trim(num2str(yStart - 1, form='(I4)')) // "-12-31 12:00:00")
! ! write data of static variables
! call var_lat%setData(wlat)
! call var_lon%setData(wlon)
! ! append data within a loop
! do i=1, ntime
! call var_time%setData(wtime(i), start=(/i/))
! call var_data%setData(wdata(:,:,i), start=(/1,1,i/))
! end do
! add some more variable attributes
call var_data%setAttribute("units", "mm/d")
call var_data%setAttribute("scaling", 0.1_dp)
......@@ -63,8 +53,6 @@ CONTAINS
! add global attributes
call nc%setAttribute("Author", trim(author_name))
end subroutine open_netcdf
end module mo_write
......@@ -58,10 +58,8 @@ end module mainVar
! RUN Control
!********************************************
module runControl
use mo_kind, only: i4, sp,dp
use mo_kind, only: i4, sp
! timer
real(sp) :: RunTime = -9.0 ! Initialization
real(sp) :: TA(2) ! Running times
character(256) :: DataPathIn !
character(256) :: DataPathOut
character(256) :: DataPathDEM
......@@ -69,12 +67,12 @@ module runControl
character(256) :: fNameSTA ! Station id and coordinates
character(256) :: fNameVARIO ! file name of variogram parameters
character(256) :: prefix ! prefix data file
character(3) :: outputformat ! nc or bin as output format
integer(i4) :: nBlocks ! number of interpolation blocks
logical :: flagEDK ! estimate EDK (T/F)
logical :: flagVario ! estimate variogran (T/F)
integer(i4) :: flagVarTyp ! type of variable 1 - peric, 2 - temp
integer(i4) :: flagMthTyp ! type of interp. method
logical :: correctNeg ! correct negative interpolated values
logical :: distZero ! values further away than dist threshold are set to zero
integer(i4) :: interMth ! interp. method
end module runControl
module kriging
......@@ -127,7 +125,7 @@ module NetCDFVar
use mo_kind, only : dp
implicit none
! Variables for NetCDF writing
character(256) :: fName_netCDF_Out ! File Name out
character(256) :: fileOut ! File Name out
real(dp), dimension(:), allocatable, target :: yCoor ! GK4 (DHDN3-zone 4) easting
real(dp), dimension(:), allocatable, target :: xCoor ! GK4 (DHDN3-zone 4) northing
real(dp), dimension(:,:), allocatable, target :: lons ! WGS84 lons
......@@ -136,4 +134,5 @@ module NetCDFVar
character(256) :: variable_name ! name of netcdf variable
character(256) :: variable_unit ! unit of netcdf variable
character(256) :: variable_long_name ! long name of netcdf variable
character(256) :: author_name ! author name of netcdf file
end module NetCDFVar
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