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

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

updated write of netcdf to mo_netcdf.f90; BUG on MacOS in the linear algebra solver; not running

parent 429c0d0f
*.g90
*/*.debug
*/*.release
\ No newline at end of file
......@@ -3,20 +3,28 @@
flagMthTyp = 1
flagVarTyp = 1
noDataValue = -9
DataPathIn = "/data/stohyd/data/processed/Germany/DWD/pre/1900-2010/"
fNameDEM = "/data/stohyd/data/processed/Germany/DEM/dkg/blocks/1.asc"
DataPathOut = "/home/zink/data/edk_test/edk/"
fNameSTA = "/data/stohyd/data/processed/Germany/DWD/pre/1900-2010/StationLuT_pre.txt"
DataPathIn = "check/pre_data/"
fNameDEM = "check/dem.asc"
DataPathOut = "output/"
FileOut = "test.nc"
fNameSTA = "check/pre_data/Stations_in_study_domain.txt"
cellFactor = 40
DataConvertFactor = 1.0e-1
flagEDK = .TRUE.
yStart = 1960
yEnd = 1960
yStart = 1990
mStart = 1
dStart = 1
yEnd = 1990
mEnd = 12
dEnd = 31
maxDist = 1.5e5
flagVario = .FALSE.
flagVario = .False.
vType = 1
nParam = 3
fNameVario = "/home/zink/fortran_code/edk/variogram_param_apriori.txt" ! path+filename
dh = 1.0e3 ! binsize for variogram
hMax = 1.5e5 ! max distance h
fNameVario = "output/varFit.txt"
dh = 1.0e3 ! binsize for variogram [m]
hMax = 1.5e5 ! max distance h [m]
! nc output specification
author_name = 'Stephan Thober'
variable_name = 'pre'
/ !END*******Main Config***********
......@@ -33,6 +33,7 @@ subroutine EDK(jd,k)
real(dp) :: sumLamda
!
! check DEM
print *, 'ha...'
if (nint(cell(k)%h) == grid%nodata_value ) then
cell(k)%z = gridMeteo%nodata_value
return
......@@ -102,11 +103,13 @@ subroutine EDK(jd,k)
B(nNmax+1) = 1.0_dp
B(nNmax+2) = cell(k)%h
print *, 'hu...'
! NOTE: only the upper triangular matrix is needed!
! call D_LSASF (A, B, X)
! call gesv(A, B) ! EVE CentOS 7
call dgesv(A, B) ! MacOS Sierra: Accelerate Framework
X = B
print *, 'hu...'
!
! The BLUE of z is then:
cell(k)%z = 0.
......
......@@ -21,8 +21,6 @@ program ED_Kriging
use mo_kind , only: i4, dp
use mo_julian , only: NDAYS, NDYIN
use mo_write_fluxes_states , only: WriteFluxStateInit, WriteFluxState, CloseFluxState_file
use runControl , only: outputformat, & ! outputformat either 'nc' or 'bin'
flagEDK, flagMthTyp, & ! flag for activate kriging, flag for 'OK' or 'EDK'
flagVarTyp, & ! pre or temp
......@@ -32,22 +30,25 @@ program ED_Kriging
nCell ! number of cells
use mo_setVario , only: setVario
use mo_netcdf , only: NcDataset, NcVariable
use mo_write , only: open_netcdf
use kriging
implicit none
character(256) :: fname
character(256) :: author_name
character(256) :: vname_data
integer(i4) :: icell ! loop varaible for cells
integer(i4) :: jDay ! loop variable - current julian day
integer(i4) :: doy ! day of year
integer(i4) :: year, month, day ! current date
integer(i4) :: netcdfid ! id of netcdf files
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
call Timer
call ReadDataMain
call ReadDataMain(fname, author_name, vname_data)
! read DEM
call ReadDEM
......@@ -66,47 +67,48 @@ program ED_Kriging
if (flagVario) call WriteDataMeteo(0,0,2)
!
if (flagEDK) then
! open netcdf if necessary
call open_netcdf(fname, ncols, nrows, nc_out, nc_data, nc_time)
timeloop: do jday = jStart, jEnd
call NDYIN(jday, day, month, year)
doy = jday - NDAYS(1,1,year) + 1
print *, 'YEAR: ',year, 'DOY: ', doy
ncellsloop: do iCell = 1, nCell
! interploation
select case (flagMthTyp)
case (1)
call EDK(jday,iCell)
case (2)
call OK(jday,iCell)
end select
end do ncellsloop
! correct precipitation values
if (flagVarTyp == 1) then
where ((cell(:)%z .LT. 0.0_sp) .AND. (cell(:)%z .GT. real(grid%nodata_value, sp)) )
cell(:)%z = 0.0_sp
end where
end if
! write output
allocate(tmp_array(gridMeteo%nrows, gridMeteo%ncols)); tmp_array=real(grid%nodata_value, dp)
tmp_array = real(reshape(cell(:)%z,(/gridMeteo%nrows, gridMeteo%ncols/)), dp)
! call WriteFluxState((jday-jStart+1), netcdfid, transpose(tmp_array))
call nc_time%setData(jday, start=(/jday/))
call nc_data%setData(tmp_array, start=(/1, 1, jday/))
deallocate(tmp_array)
end do timeloop
! close netcdf if necessary
call nc_out%close()
! open netcdf if necessary
call open_netcdf(fname, author_name, vname_data, grid%ncols, grid%nrows, yStart, nc_out, nc_data, nc_time)
timeloop: do jday = jStart, jEnd
call NDYIN(jday, day, month, year)
doy = jday - NDAYS(1,1,year) + 1
print *, 'YEAR: ',year, 'DOY: ', doy
ncellsloop: do iCell = 1, nCell
! interploation
select case (flagMthTyp)
case (1)
print *, 'ha...'
call EDK(jday,iCell)
case (2)
call OK(jday,iCell)
end select
end do ncellsloop
print *, 'ha...'
! correct precipitation values
if (flagVarTyp == 1) then
where ((cell(:)%z .LT. 0.0_sp) .AND. (cell(:)%z .GT. real(grid%nodata_value, sp)) )
cell(:)%z = 0.0_sp
end where
end if
! write output
allocate(tmp_array(gridMeteo%nrows, gridMeteo%ncols)); tmp_array=real(grid%nodata_value, dp)
tmp_array = real(reshape(cell(:)%z,(/gridMeteo%nrows, gridMeteo%ncols/)), dp)
call nc_time%setData(jday, start=(/jday/))
call nc_data%setData(tmp_array, start=(/1, 1, jday/))
deallocate(tmp_array)
end do timeloop
! close netcdf if necessary
call nc_out%close()
end if
! deallocate memory
......
......@@ -134,7 +134,7 @@ subroutine EmpVar(jd, flagMax)
if (Nh(k) > 0) then
ni=ni+1
! Classsical
gamma(k,2)=gamma(k,2)/2.0_dp/real(Nh(k), dp)
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
......@@ -149,7 +149,7 @@ subroutine EmpVar(jd, flagMax)
!
! scaling
if (flagMax) gmax=maxval(gamma, dim=1)
do k=1,nbins
if (gamma(k,1) > 0.0_dp) then
gamma(k,1) = gamma(k,1)/gmax(1)
......
......@@ -7,7 +7,7 @@
! Created Sa 21.03.2006
! Last Update Sa 11.06.2010
!**************************************************************************
subroutine ReadDataMain
subroutine ReadDataMain(FileOut, author_name, variable_name)
use mainVar
use mo_kind, only: i4, dp
use mo_julian , only: NDAYS
......@@ -15,22 +15,25 @@ subroutine ReadDataMain
use kriging
use VarFit
use runControl
use NetCDFVar, only: variable_name, variable_unit, variable_long_name
use NetCDFVar, only: variable_unit, variable_long_name
implicit none
!
integer(i4) :: i, ios
character(256) :: dummy, fileName
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
!
!===============================================================
! namelist definition
!===============================================================
!
namelist/mainVars/flagMthTyp, flagVarTyp, noDataValue, DataPathIn, fNameDEM, &
DataPathOut, fNameSTA, cellFactor, outputformat, DataConvertFactor, flagEDK, &
variable_name, variable_unit, variable_long_name, &
DataPathOut, FileOut, fNameSTA, cellFactor, outputformat, DataConvertFactor, flagEDK, &
author_name, variable_name, variable_unit, variable_long_name, &
yStart, mStart, dStart, yEnd, mEnd, dEnd, maxDist, flagVario, vType, nParam, &
fNameVario, dh, hMax
fNameVario, dh, hMax, author_name, variable_name
!
! -----------------------------------------------------------------------
! MAIN.DAT
......@@ -38,6 +41,7 @@ subroutine ReadDataMain
open(unit=10, file='main.dat', STATUS='OLD', ACTION='read')
read(10, nml=mainVars)
close(10)
FileOut = trim(DataPathOut) // trim(FileOut)
!
! *************************************
! Read initial control parameters
......
This diff is collapsed.
!> \file mo_set_netcdf_outputs.f90
!> \brief Defines the structure of the netCDF to write the output in.
!> \details All output variables are initialized for the NetCDF.
!> \authors Matthias Zink
!> \date Apr 2013
MODULE mo_set_netcdf_outputs
! This module defines the NetCDF variables.
! Written Matthias Zink, Apr 2013
USE mo_kind, ONLY: i4
IMPLICIT NONE
PRIVATE
PUBLIC :: set_netCDF
! ------------------------------------------------------------------
CONTAINS
! ------------------------------------------------------------------
! NAME
! WriteFluxStateInit
!> \brief Runs mhm with a specific parameter set and returns required variables, e.g. runoff.
!> \details Runs mhm with a specific parameter set and returns required variables, e.g. runoff.
! INTENT(IN)
! None
! INTENT(INOUT)
! None
! INTENT(OUT)
! None
! INTENT(IN), OPTIONAL
! None
! INTENT(INOUT), OPTIONAL
! None
! INTENT(OUT), OPTIONAL
! None
! RETURN
! None
! RESTRICTIONS
!
! EXAMPLE
!
! LITERATURE
! HISTORY
!> \author Matthias Zink
!> \date Apr 2013
! Modified,
subroutine set_netCDF(NoNetcdfVars, nrows, ncols)
!
use mo_NCWrite , only: V, Dnc, ndims, nVars
use netcdf , only: NF90_UNLIMITED, NF90_CHAR, NF90_INT, NF90_DOUBLE
implicit none
!
integer(i4), intent(in) :: NoNetcdfVars
integer(i4), intent(in) :: nrows
integer(i4), intent(in) :: ncols
!
integer(i4) :: i
!
! define output file
!
! define parameters
nDims = 3 ! nr. dim types
nVars = 5 + NoNetcdfVars ! total nr. var to print
!
! allocate arrays
if (allocated(Dnc) ) deallocate ( Dnc)
if (allocated(V) ) deallocate ( V )
allocate ( Dnc(nDims) )
allocate ( V(nVars) )
!
! define dimensions (coordinates) => corresponding variable must be created
i = 1
Dnc(i)%name = "easting"
Dnc(i)%len = nrows
!
i = 2
Dnc(i)%name = "northing"
Dnc(i)%len = ncols
!
i = 3
Dnc(i)%name = "time"
Dnc(i)%len = NF90_UNLIMITED
!
!
! DIMENSION VARIABLES
i = 1
V(i)%name = Dnc(i)%name
V(i)%xType = NF90_DOUBLE
V(i)%nLvls = 0
V(i)%nSubs = 0
V(i)%nDims = 1
V(i)%dimTypes = (/i,0,0,0,0/)
V(i)%wFlag = .true.
! attributes (other possibilities: add_offset, valid_min, valid_max)
V(i)%nAtt = 2
!
V(i)%att(1)%name = "units"
V(i)%att(1)%xType = NF90_CHAR
V(i)%att(1)%nValues= 1
V(i)%att(1)%values = "m"
!
V(i)%att(2)%name = "long_name"
V(i)%att(2)%xType = NF90_CHAR
V(i)%att(2)%nValues= 1
V(i)%att(2)%values = "x-coordinate in cartesian coordinates GK4"
!
i = 2
V(i)%name = Dnc(i)%name
V(i)%xType = NF90_DOUBLE
V(i)%nLvls = 0
V(i)%nSubs = 0
V(i)%nDims = 1
V(i)%dimTypes = (/i,0,0,0,0/)
! printing
V(i)%wFlag = .true.
! attributes
V(i)%nAtt = 2
!
V(i)%att(1)%name = "units"
V(i)%att(1)%xType = NF90_CHAR
V(i)%att(1)%nValues= 1
V(i)%att(1)%values = "m"
!
V(i)%att(2)%name = "long_name"
V(i)%att(2)%xType = NF90_CHAR
V(i)%att(2)%nValues= 1
V(i)%att(2)%values = "y-coordinate in cartesian coordinates GK4"
!
i = 3
V(i)%name = Dnc(i)%name
V(i)%xType = NF90_INT
V(i)%nLvls = 0
V(i)%nSubs = 0
V(i)%nDims = 1
V(i)%dimTypes = (/i,0,0,0,0/)
V(i)%wFlag = .true.
! attributes
V(i)%nAtt = 2
!
V(i)%att(1)%name = "units"
V(i)%att(1)%xType = NF90_CHAR
V(i)%att(1)%nValues= 1
!
V(i)%att(2)%name = "long_name"
V(i)%att(2)%xType = NF90_CHAR
V(i)%att(2)%nValues= 1
V(i)%att(2)%values = "time"
!
i = 4
V(i)%name = "lon"
V(i)%xType = NF90_DOUBLE
V(i)%nLvls = 1
V(i)%nSubs = 1
V(i)%nDims = 2
V(i)%dimTypes = (/1,2,0,0,0/)
V(i)%wFlag = .true.
! attributes
V(i)%nAtt = 2
!
V(i)%att(1)%name = "units"
V(i)%att(1)%xType = NF90_CHAR
V(i)%att(1)%nValues= 1
V(i)%att(1)%values = "degrees_east"
!
V(i)%att(2)%name = "long_name"
V(i)%att(2)%xType = NF90_CHAR
V(i)%att(2)%nValues= 1
V(i)%att(2)%values = "longitude"
!
i = 5
V(i)%name = "lat"
V(i)%xType = NF90_DOUBLE
V(i)%nLvls = 1
V(i)%nSubs = 1
V(i)%nDims = 2
V(i)%dimTypes = (/1,2,0,0,0/)
V(i)%wFlag = .true.
! attributes
V(i)%nAtt = 2
!
V(i)%att(1)%name = "units"
V(i)%att(1)%xType = NF90_CHAR
V(i)%att(1)%nValues= 1
V(i)%att(1)%values = "degrees_north"
!
V(i)%att(2)%name = "long_name"
V(i)%att(2)%xType = NF90_CHAR
V(i)%att(2)%nValues= 1
V(i)%att(2)%values = "latitude"
!
! ******************************
! dynamic variables data
! ******************************
do i = 6, NoNetcdfVars + 5
V(i)%xType = NF90_DOUBLE
V(i)%nLvls = 1
V(i)%nSubs = 1
V(i)%nDims = 3
V(i)%dimTypes = (/1,2,3,0,0/)
! printing
V(i)%wFlag = .true.
! pointer
! during running time
! attributes
V(i)%nAtt = 6
!
V(i)%att(1)%name = "units"
V(i)%att(1)%xType = NF90_CHAR
V(i)%att(1)%nValues= 1
!
V(i)%att(2)%name = "long_name"
V(i)%att(2)%xType = NF90_CHAR
V(i)%att(2)%nValues= 1
!
V(i)%att(3)%name = "scale_factor"
V(i)%att(3)%xType = NF90_DOUBLE
V(i)%att(3)%nValues= 1
!
V(i)%att(4)%name = "_FillValue"
V(i)%att(4)%xType = NF90_DOUBLE
V(i)%att(4)%nValues= 1
!
V(i)%att(5)%name = "missing_value"
V(i)%att(5)%xType = NF90_DOUBLE
V(i)%att(5)%nValues= 1
! !
V(i)%att(6)%name = "coordinates"
V(i)%att(6)%xType = NF90_CHAR
V(i)%att(6)%nValues= 1
end do
!
!
end subroutine set_netCDF
END MODULE mo_set_netcdf_outputs
......@@ -8,16 +8,20 @@ module mo_write
CONTAINS
subroutine open_netcdf(fname, ncols, nrows, nc_id, var_id, var_time)
subroutine open_netcdf(fname, author_name, vname_data, ncols, nrows, year_start, 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
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
......@@ -33,13 +37,13 @@ CONTAINS
dim_time = nc%setDimension("time", -1)
! create variables
var_time = nc%setVariable(vname_time, "i32", (/dim_time/))
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_x, dim_y, dim_time/))
! add some variable attributes
call var_time%setAttribute("units", "days since " // trim(num2str(yearstart - 1, form='(I4)')) // "-12-31 12:00:00")
call var_time%setAttribute("units", "days since " // trim(num2str(year_start - 1, form='(I4)')) // "-12-31 12:00:00")
! ! write data of static variables
! call var_lat%setData(wlat)
......
!> \file mo_write_fluxes_states.f90
!> \brief Creates NetCDF output for different fluxes and state variabels of mHM.
!> \details NetCDF is first initialized and later on variables are put to the NetCDF.
!> \authors Matthias Zink
!> \date Apr 2013
MODULE mo_write_fluxes_states
! This module creates the output for mHM.
! Written Matthias Zink, Apr 2012
USE mo_kind, ONLY: i4, dp
IMPLICIT NONE
PRIVATE
PUBLIC :: WriteFluxStateInit
PUBLIC :: WriteFluxState
PUBLIC :: CloseFluxState_file
! ------------------------------------------------------------------