Browse Source

Added bias correction in flux_calculator

This feature uses python calling from Fortran
1.03.00
Sven Karsten 5 months ago
parent
commit
57119cd371
  1. 11
      src/bias_corrections.F90
  2. 1
      src/flux_calculator.F90
  3. 3
      src/flux_calculator_basic.F90
  4. 20
      src/flux_calculator_calculate.F90
  5. 78
      src/pyfort/call_python.f90
  6. 13
      src/pyfort/datetime_helpers.py
  7. 111
      src/pyfort/pyfort_src.py

11
src/bias_corrections.F90

@ -19,7 +19,7 @@ ENUM, BIND(C)
ENDENUM
! initializes names for corrextion -> corresponds to out variable that is corrected
CHARACTER(len=8), PARAMETER, DIMENSION(E_N_CORRECTIONS) :: corrections_names = [ &
CHARACTER(len=10), PARAMETER, DIMENSION(E_N_CORRECTIONS) :: corrections_names = [ &
'mass_evap' &
]
@ -190,6 +190,9 @@ SUBROUTINE initialize_bias_corrections(filename, grid_offset, grid_size)
! 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
@ -202,9 +205,9 @@ SUBROUTINE initialize_bias_corrections(filename, grid_offset, grid_size)
correction_filename = 'corrections/'//TRIM(corrections_names(i))//'-'//TRIM(yj)//'.nc'
istatus = nf90_open(TRIM(filename), NF90_NOWRITE, ncfileid)
istatus = nf90_open(TRIM(correction_filename), NF90_NOWRITE, ncfileid)
IF (istatus /= NF90_NOERR) THEN
WRITE(*,*) 'Could not open ', TRIM(filename), ' for bias correction. Unset correction.'
WRITE(*,*) 'Could not open ', TRIM(correction_filename), ' for bias correction. Unset correction.'
CYCLE
ENDIF
@ -232,7 +235,7 @@ SUBROUTINE initialize_bias_corrections(filename, grid_offset, grid_size)
istatus = nf90_close(ncfileid)
IF (istatus /= NF90_NOERR) THEN
WRITE(*,*) 'Could not close ', TRIM(filename), 'for bias correction.'
WRITE(*,*) 'Could not close ', TRIM(correction_filename), 'for bias correction.'
CYCLE
ENDIF

1
src/flux_calculator.F90

@ -858,6 +858,7 @@ ENDIF
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