Browse Source

Merge branch 'experiments/bias-correction'

1.03.00
Sven Karsten 5 months ago
parent
commit
4e8d0810e1
  1. 11
      build_hlrnb.sh
  2. 11
      build_hlrng.sh
  3. 257
      src/bias_corrections.F90
  4. 9
      src/flux_calculator.F90
  5. 3
      src/flux_calculator_basic.F90
  6. 20
      src/flux_calculator_calculate.F90
  7. 78
      src/pyfort/call_python.f90
  8. 13
      src/pyfort/datetime_helpers.py
  9. 111
      src/pyfort/pyfort_src.py

11
build_hlrnb.sh

@ -8,6 +8,8 @@ module unload impi
module load intel/18.0.6
module load impi/2018.5
module load anaconda3/2019.03
FC=mpiifort
AR=ar
IOW_ESM_ROOT="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )/../.."
@ -59,12 +61,19 @@ $FC -c $FFLAGS ../src/flux_lib/*.F90
rm flux_library.a
$AR rv flux_library.a *.o
cp -r ${PWD}/../src/pyfort/* .
python3 pyfort_src.py "$PWD"
$FC -c -Wl,-rpath,"${PWD}" -L"${PWD}" -lpyfort call_python.f90
$FC -c $FFLAGS ../src/bias_corrections.F90 -I${IOW_ESM_NETCDF_INCLUDE} $LIBS
$FC -c $FFLAGS ../src/flux_calculator_basic.F90
$FC -c $FFLAGS ../src/flux_calculator_prepare.F90
$FC -c $FFLAGS ../src/flux_calculator_calculate.F90
$FC -c $FFLAGS ../src/flux_calculator_parse_arg.F90
$FC -c $FFLAGS ../src/flux_calculator_io.F90 -I${IOW_ESM_NETCDF_INCLUDE} $LIBS
$FC -c $FFLAGS ../src/flux_calculator_create_namcouple.F90
$FC $FFLAGS -o ../"${bin_dir}"/flux_calculator ../src/flux_calculator.F90 flux_calculator_basic.o flux_calculator_prepare.o flux_calculator_calculate.o flux_calculator_io.o flux_calculator_parse_arg.o flux_calculator_create_namcouple.o flux_library.a $INCLUDES $LIBS -Wl,-rpath,${IOW_ESM_NETCDF_LIBRARY}
$FC $FFLAGS -o ../"${bin_dir}"/flux_calculator ../src/flux_calculator.F90 flux_calculator_basic.o flux_calculator_prepare.o flux_calculator_calculate.o flux_calculator_io.o flux_calculator_parse_arg.o flux_calculator_create_namcouple.o bias_corrections.o call_python.o flux_library.a \
$INCLUDES $LIBS -Wl,-rpath,"$PWD",-rpath,${IOW_ESM_NETCDF_LIBRARY},-rpath=/sw/tools/python/anaconda3/2019.03/skl/lib -L/sw/tools/python/anaconda3/2019.03/skl/lib -lpython3.7m -L"$PWD" -lpyfort
cd -

11
build_hlrng.sh

@ -8,6 +8,8 @@ module unload impi
module load intel/18.0.5
module load impi/2018.5
module load anaconda3/2019.03
FC=mpiifort
AR=ar
IOW_ESM_ROOT="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )/../.."
@ -59,12 +61,19 @@ $FC -c $FFLAGS ../src/flux_lib/*.F90
rm flux_library.a
$AR rv flux_library.a *.o
cp -r ${PWD}/../src/pyfort/* .
python3 pyfort_src.py "$PWD"
$FC -c -Wl,-rpath,"${PWD}" -L"${PWD}" -lpyfort call_python.f90
$FC -c $FFLAGS ../src/bias_corrections.F90 -I${IOW_ESM_NETCDF_INCLUDE} $LIBS
$FC -c $FFLAGS ../src/flux_calculator_basic.F90
$FC -c $FFLAGS ../src/flux_calculator_prepare.F90
$FC -c $FFLAGS ../src/flux_calculator_calculate.F90
$FC -c $FFLAGS ../src/flux_calculator_parse_arg.F90
$FC -c $FFLAGS ../src/flux_calculator_io.F90 -I${IOW_ESM_NETCDF_INCLUDE} $LIBS
$FC -c $FFLAGS ../src/flux_calculator_create_namcouple.F90
$FC $FFLAGS -o ../"${bin_dir}"/flux_calculator ../src/flux_calculator.F90 flux_calculator_basic.o flux_calculator_prepare.o flux_calculator_calculate.o flux_calculator_io.o flux_calculator_parse_arg.o flux_calculator_create_namcouple.o flux_library.a $INCLUDES $LIBS -Wl,-rpath,${IOW_ESM_NETCDF_LIBRARY}
$FC $FFLAGS -o ../"${bin_dir}"/flux_calculator ../src/flux_calculator.F90 flux_calculator_basic.o flux_calculator_prepare.o flux_calculator_calculate.o flux_calculator_io.o flux_calculator_parse_arg.o flux_calculator_create_namcouple.o bias_corrections.o call_python.o flux_library.a \
$INCLUDES $LIBS -Wl,-rpath,"$PWD",-rpath,${IOW_ESM_NETCDF_LIBRARY},-rpath=/sw/tools/python/anaconda3/2019.03/skl/lib -L/sw/tools/python/anaconda3/2019.03/skl/lib -lpython3.7m -L"$PWD" -lpyfort
cd -

257
src/bias_corrections.F90

@ -0,0 +1,257 @@
MODULE bias_corrections
USE netcdf, ONLY : &
nf90_open, &
nf90_close, &
nf90_get_var, &
nf90_inq_varid, &
nf90_get_att, &
NF90_NOWRITE, &
NF90_NOERR
IMPLICIT NONE
PUBLIC initialize_bias_corrections
ENUM, BIND(C)
ENUMERATOR :: E_MASS_EVAP_CORRECTION = 1
ENUMERATOR :: E_N_CORRECTIONS = 1
ENDENUM
! initializes names for corrextion -> corresponds to out variable that is corrected
CHARACTER(len=10), PARAMETER, DIMENSION(E_N_CORRECTIONS) :: corrections_names = [ &
'mass_evap' &
]
INTEGER :: &
init_date
REAL, ALLOCATABLE :: &
corrections(:,:,:) ! (variable, month, space)
LOGICAL :: &
lcorrections(E_N_CORRECTIONS) = .FALSE.
CONTAINS
SUBROUTINE process_input_corrections (n, errstat)
! Parameter list:
INTEGER, INTENT (IN) :: &
n ! Unit number for Namelist INPUT file
INTEGER, INTENT (OUT) :: &
errstat ! error status variable
! Local variables:
INTEGER :: &
i
! Define the namelist group
NAMELIST /correctionsctl/ init_date, lcorrections
!------------------------------------------------------------------------------
!- End of header -
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
!- Begin SUBROUTINE input_oasisctl
!------------------------------------------------------------------------------
errstat = 0
!------------------------------------------------------------------------------
!- Section 3: Input of the namelist values
!------------------------------------------------------------------------------
READ (n, correctionsctl, IOSTAT=errstat)
IF (errstat /= 0) THEN
errstat = 1
DO i = 1, E_N_CORRECTIONS
lcorrections(i) = .FALSE.
ENDDO
WRITE(*,*) "Could not read INPUT_BIAS set lcorrections to default: ", lcorrections
RETURN
ENDIF
!------------------------------------------------------------------------------
!- Section 4: Check values for errors and consistency
!------------------------------------------------------------------------------
DO i = 1, E_N_CORRECTIONS
IF ( lcorrections(i) /= .FALSE. .and. lcorrections(i) /= .TRUE.) THEN
WRITE (*,*) ' ERROR *** Wrong value ', lcorrections(i) , ' for correction ', i
errstat = 2
WRITE (*,*) ' ERROR *** Error while checking values of the namelist oasisctl *** '
RETURN
ENDIF
ENDDO
!------------------------------------------------------------------------------
!- Section 6: Output of the namelist variables and their default values
!------------------------------------------------------------------------------
WRITE (*,*) 'Apply bias corrections for:'
DO i = 1, E_N_CORRECTIONS
WRITE (*,*) corrections_names(i), lcorrections(i)
ENDDO
!------------------------------------------------------------------------------
!- End of the Subroutine
!------------------------------------------------------------------------------
END SUBROUTINE process_input_corrections
SUBROUTINE read_namelist_corrections(filename, errstat)
CHARACTER (LEN= *), INTENT(IN) :: &
filename ! error message
INTEGER, INTENT(OUT) :: errstat
INTEGER :: n
n = 10
! -----------------------------------------------------------------
! 1 Open NAMELIST-INPUT file
! ----------------------------------------------------------------
errstat = 0
WRITE (*,*) ' INPUT OF THE NAMELIST FOR BIAS CORRECTIONS'
OPEN (n, FILE=filename, FORM='FORMATTED', STATUS='UNKNOWN', IOSTAT=errstat)
IF (errstat /= 0) THEN
errstat = 1
WRITE(*,*) ' ERROR *** Error while opening file '//filename//' *** '
RETURN
ENDIF
! -----------------------------------------------------------------
! 2 read the NAMELIST group oasisctl
! ----------------------------------------------------------------
CALL process_input_corrections (n, errstat)
IF (errstat /= 0) THEN
WRITE (*,*) ' ERROR *** Wrong values occured in NAMELIST group /correctionsctl/ *** '
errstat = 2
RETURN
ENDIF
! -----------------------------------------------------------------
! 3 Close NAMELIST-INPUT file
! ----------------------------------------------------------------
CLOSE (n, STATUS='KEEP', IOSTAT=errstat)
IF (errstat /= 0) THEN
WRITE(*,*) ' ERROR *** while closing file '//filename//'*** '
errstat = 4
ENDIF
!------------------------------------------------------------------------------
!- End of the Subroutine
!------------------------------------------------------------------------------
END SUBROUTINE read_namelist_corrections
SUBROUTINE initialize_bias_corrections(filename, grid_offset, grid_size)
CHARACTER (LEN= *), INTENT(IN) :: &
filename ! name of input file
INTEGER, INTENT(IN) :: &
grid_offset, &
grid_size
INTEGER :: &
istatus, & ! NetCDF status
ncfileid, ncvarid, & ! NetCDF IDs
nerror, &
i, j
REAL:: &
fillvalue
CHARACTER(LEN=128) :: correction_filename
CHARACTER(LEN=2) :: yj
! Read namelist for corrections
CALL read_namelist_corrections(filename, nerror)
! allocate resources for each process
ALLOCATE(corrections(E_N_CORRECTIONS, 12, grid_size))
! intialize with zero
corrections(:,:,:) = 0.0
! Read in correction for each month
DO i = 1, E_N_CORRECTIONS
IF (.NOT. lcorrections(i)) THEN
CYCLE
ENDIF
DO j = 1, 12
WRITE (yj,'(I2.2)') j
yj = ADJUSTL(yj)
correction_filename = 'corrections/'//TRIM(corrections_names(i))//'-'//TRIM(yj)//'.nc'
istatus = nf90_open(TRIM(correction_filename), NF90_NOWRITE, ncfileid)
IF (istatus /= NF90_NOERR) THEN
WRITE(*,*) 'Could not open ', TRIM(correction_filename), ' for bias correction. Unset correction.'
CYCLE
ENDIF
istatus = nf90_inq_varid(ncfileid, TRIM(corrections_names(i)) , ncvarid)
IF (istatus /= NF90_NOERR) THEN
WRITE(*,*) 'Could not get varid for variable '//TRIM(corrections_names(i))//'. Unset correction.'
CYCLE
ENDIF
istatus = nf90_get_var(ncfileid, ncvarid, corrections(i,j,:), &
(/ grid_offset/), &
(/ grid_size/))
IF (istatus /= NF90_NOERR) THEN
WRITE(*,*) 'Could not get variable '//TRIM(corrections_names(i))//'. Unset correction.'
corrections(i,j,:) = 0.0
CYCLE
ENDIF
istatus = nf90_get_att(ncfileid, ncvarid, "_FillValue", fillvalue)
IF (istatus /= NF90_NOERR) THEN
WRITE(*,*) 'Could not get fill value. Unset correction.'
corrections(i,j,:) = 0.0
CYCLE
ENDIF
istatus = nf90_close(ncfileid)
IF (istatus /= NF90_NOERR) THEN
WRITE(*,*) 'Could not close ', TRIM(correction_filename), 'for bias correction.'
CYCLE
ENDIF
WHERE (corrections(i,j,:) == fillvalue) corrections(i,j,:) = 0.0
ENDDO
ENDDO
WRITE (*,*) "Read in bias correction fields."
END SUBROUTINE initialize_bias_corrections
SUBROUTINE finalize_bias_corrections
DEALLOCATE(corrections)
END SUBROUTINE finalize_bias_corrections
END MODULE bias_corrections

9
src/flux_calculator.F90

@ -15,6 +15,8 @@ PROGRAM flux_calculator
USE flux_calculator_parse_arg
USE flux_calculator_create_namcouple
USE bias_corrections, ONLY : initialize_bias_corrections
! Use OASIS communication library
USE mod_oasis
@ -175,6 +177,8 @@ PROGRAM flux_calculator
WRITE(*,nml=input)
CLOSE(10)
! Initialize the idx_???? variables which store the index of a variable name
CALL init_varname_idx
@ -844,12 +848,17 @@ ENDIF
ENDIF
WRITE (w_unit,*) 'Finished initialization.'
!### if available, initialize bias corrections
CALL initialize_bias_corrections('flux_calculator.nml', grid_offset(1), grid_size(1))
!###############################################################################
!# STEP 2: TIME LOOP #
!###############################################################################
DO n_timestep = 1,num_timesteps
current_time = (n_timestep - 1) * timestep
current_step_time = current_time ! update global time counter
IF (verbosity_level >= VERBOSITY_LEVEL_STANDARD) THEN
WRITE (w_unit,*) 'Time since start = ',current_time,' seconds.'
CALL FLUSH(w_unit)

3
src/flux_calculator_basic.F90

@ -121,6 +121,9 @@ MODULE flux_calculator_basic
TYPE(realarray) :: weight
END TYPE sparse_regridding_matrix
! time of current step in seconds from the start of this instance (updated main loop)
INTEGER :: current_step_time = 0
CONTAINS
SUBROUTINE add_input_field(myname, my_letter, surface_type, which_grid, local_field, num_input_fields, input_field)

20
src/flux_calculator_calculate.F90

@ -3,6 +3,10 @@ MODULE flux_calculator_calculate
use flux_calculator_basic
use flux_library
use bias_corrections, only: lcorrections, corrections, E_MASS_EVAP_CORRECTION, init_date
use, intrinsic :: iso_c_binding
use call_python, only : get, set, call_function
IMPLICIT NONE
@ -57,6 +61,16 @@ MODULE flux_calculator_calculate
CHARACTER(len=4), PARAMETER :: myvarname = 'MEVA'
CHARACTER(len=20) :: method
INTEGER :: i, j
INTEGER :: current_month(1)
! if we want to apply a correction we need to know the month (use python interface for that)
IF (lcorrections(E_MASS_EVAP_CORRECTION)) THEN
CALL set("init_date", init_date)
CALL set("seconds", current_step_time)
call call_function("datetime_helpers", "get_current_date")
call get("current_month", current_month)
!WRITE (*,*) "Current month = ", current_month(1)
ENDIF
DO i=1,num_surface_types
method = methods(my_bottom_model,i)
@ -94,6 +108,12 @@ MODULE flux_calculator_calculate
local_field(i,1)%var(idx_VATM)%field(j))
ENDDO
ENDIF
IF (lcorrections(E_MASS_EVAP_CORRECTION)) THEN
DO j=1,grid_size(1)
local_field(i,1)%var(idx_MEVA)%field(j) = local_field(i,1)%var(idx_MEVA)%field(j) + corrections(E_MASS_EVAP_CORRECTION, current_month(1), j)
ENDDO
ENDIF
ENDIF
ENDDO
!CALL average_across_surface_types(1,idx_MEVA,num_surface_types,grid_size,local_field)

78
src/pyfort/call_python.f90

@ -0,0 +1,78 @@
module call_python
use, intrinsic :: iso_c_binding
implicit none
private
public :: get
public :: set
public :: call_function
contains
subroutine check(ret)
integer :: ret
if (ret /= 0) stop -1
end subroutine check
subroutine set(variable_name, variable_value)
interface
function set_(variable_name_c, variable_value_c) result(y) bind (c, name='set')
use iso_c_binding
character(c_char) :: variable_name_c
integer(c_int) :: variable_value_c
integer(c_int) :: y
end function set_
end interface
character(len=*) :: variable_name
integer(c_int) :: variable_value
character(kind=c_char, len=128) :: variable_name_c
integer(c_int) :: variable_value_c
variable_name_c = trim(variable_name)//char(0)
variable_value_c = variable_value
call check(set_(variable_name_c, variable_value_c))
end subroutine set
subroutine get(variable_name, variable_value)
interface
function get_(variable_name_c, variable_value_c) result(y) bind (c, name='get')
use iso_c_binding
character(c_char) :: variable_name_c
integer(c_int) :: variable_value_c(1)
integer(c_int) :: y
end function get_
end interface
character(len=*) :: variable_name
integer(c_int) :: variable_value(:)
character(kind=c_char, len=128) :: variable_name_c
integer(c_int) :: variable_value_c(size(variable_value, 1))
variable_name_c = trim(variable_name)//char(0)
variable_value_c = variable_value
call check(get_(variable_name_c, variable_value_c))
variable_value = variable_value_c
end subroutine get
subroutine call_function(module_name, function_name)
interface
function call_function_(mod_name_c, fun_name_c) result(y) bind(c, name='call_function')
use iso_c_binding
character(c_char) mod_name_c, fun_name_c
integer(c_int) :: y
end function call_function_
end interface
character(len=*) :: module_name, function_name
character(kind=c_char, len=256) :: mod_name_c, fun_name_c
mod_name_c = trim(module_name)//char(0)
fun_name_c = trim(function_name)//char(0)
call check(call_function_(mod_name_c, fun_name_c))
end subroutine call_function
end module call_python

13
src/pyfort/datetime_helpers.py

@ -0,0 +1,13 @@
from datetime import datetime
from datetime import timedelta
def get_current_date(state):
init_date = str(*state["init_date"])
seconds = int(*state["seconds"])
dt = datetime.strptime(init_date, "%Y%m%d")
current_date = dt + timedelta(seconds=seconds)
state["current_date"] = int(current_date.strftime("%Y%m%d"))
state["current_month"] = int(current_date.strftime("%m"))
state["current_day"] = int(current_date.strftime("%d"))
state["current_year"] = int(current_date.strftime("%Y"))

111
src/pyfort/pyfort_src.py

@ -0,0 +1,111 @@
import sys
import os
try:
module_dir = str(sys.argv[1])
except:
module_dir = os.getcwd()
import cffi
ffibuilder = cffi.FFI()
header = """
extern int set(char *, int32_t *);
extern int get(char *, int32_t *);
extern int call_function(char *, char *);
"""
module = f"""
from pyfort import ffi
import numpy as np
import importlib
import os
import sys
# Create the dictionary mapping ctypes to np dtypes.
ctype2dtype = {{}}
# Integer types
for prefix in ('int', 'uint'):
for log_bytes in range(4):
ctype = '%s%d_t' % (prefix, 8 * (2**log_bytes))
dtype = '%s%d' % (prefix[0], 2**log_bytes)
#print( ctype )
#print( dtype )
ctype2dtype[ctype] = np.dtype(dtype)
# Floating point types
ctype2dtype['float'] = np.dtype('f4')
ctype2dtype['double'] = np.dtype('f8')
def asarray(ffi, ptr, shape, **kwargs):
length = np.prod(shape)
# Get the canonical C type of the elements of ptr as a string.
T = ffi.getctype(ffi.typeof(ptr).item)
#print( T )
#print( ffi.sizeof( T ) )
if T not in ctype2dtype:
raise RuntimeError("Cannot create an array for element type: %s" % T)
a = np.frombuffer(ffi.buffer(ptr, length * ffi.sizeof(T)), ctype2dtype[T]).reshape(shape, **kwargs)
return a
STATE = {{}}
@ffi.def_extern(error=1)
def get(key, c_ptr, shape=(1,)):
\"\"\"Copy the numpy array stored in STATE[key] to a pointer\"\"\"
key = ffi.string(key).decode("UTF-8")
# wrap pointer in numpy array
fortran_arr = asarray(ffi, c_ptr, shape)
# update the numpy array in place
fortran_arr[:] = STATE[key]
return 0
@ffi.def_extern(error=1)
def set(key, c_ptr, shape=(1,)):
\"\"\"Call python\"\"\"
key = ffi.string(key).decode("UTF-8")
STATE[key] = asarray(ffi, c_ptr, shape).copy()
return 0
@ffi.def_extern(error=1)
def call_function(module_name, function_name):
#pwd = os.getcwd()
#print(pwd)
#sys.path.append(pwd)
sys.path.append("{module_dir}")
module_name = ffi.string(module_name).decode("UTF-8")
function_name = ffi.string(function_name).decode("UTF-8")
mod = importlib.import_module(module_name)
# the function we want to call
fun = getattr(mod, function_name)
# call the function
# this function can edit STATE inplace
fun(STATE)
return 0
"""
with open("pyfort.h", "w") as f:
f.write(header)
ffibuilder.embedding_api(header)
ffibuilder.set_source("pyfort", r'''
#include "pyfort.h"
''',)
ffibuilder.embedding_init_code(module)
ffibuilder.emit_c_code("pyfort.c")
ffibuilder.compile(target="libpyfort.so", verbose=True)
Loading…
Cancel
Save