Commit 866d2c1d authored by Jorn Bruggeman's avatar Jorn Bruggeman
Browse files

python driver: added do_surface, do_bottom flags to get_rates

parent fa2292f0
......@@ -108,7 +108,7 @@ fabm.check_ready.argtypes = []
fabm.check_ready.restype = None
# Routine for retrieving source-sink terms for the interior domain.
fabm.get_rates.argtypes = [numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags='CONTIGUOUS')]
fabm.get_rates.argtypes = [numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags='CONTIGUOUS'), ctypes.c_int, ctypes.c_int]
fabm.get_rates.restype = None
# Routine for getting git repository version information.
......@@ -451,12 +451,12 @@ class Model(object):
self.bulk_state_variables = self.interior_state_variables
self.bulk_diagnostic_variables = self.interior_diagnostic_variables
def getRates(self):
def getRates(self, surface=True, bottom=True):
"""Returns the local rate of change in state variables,
given the current state and environment.
"""
localrates = numpy.empty_like(self.state)
fabm.get_rates(localrates)
fabm.get_rates(localrates, surface, bottom)
return localrates
def getJacobian(self,pert=None):
......
......@@ -386,27 +386,29 @@
call fabm_link_bottom_state_data(model,index,value)
end subroutine link_bottom_state_data
subroutine get_rates(pelagic_rates_) bind(c)
subroutine get_rates(pelagic_rates_, do_surface, do_bottom) bind(c)
!DIR$ ATTRIBUTES DLLEXPORT :: get_rates
real(c_double),target,intent(in) :: pelagic_rates_(*)
integer(c_int),value, intent(in) :: do_surface, do_bottom
real(c_double),pointer :: pelagic_rates(:)
real(rk) :: ext
if (.not.associated(column_depth)) call driver%fatal_error('get_rates', &
'Value for environmental dependency '//trim(environment_names(index_column_depth))// &
' must be provided before calling get_rates.')
call fabm_get_light_extinction(model,ext)
call fabm_get_light(model)
call c_f_pointer(c_loc(pelagic_rates_),pelagic_rates, &
(/size(model%state_variables)+size(model%surface_state_variables)+size(model%bottom_state_variables)/))
pelagic_rates = 0.0_rk
call fabm_do_surface(model,pelagic_rates(1:size(model%state_variables)), &
if (int2logical(do_surface)) call fabm_do_surface(model,pelagic_rates(1:size(model%state_variables)), &
pelagic_rates(size(model%state_variables)+1:size(model%state_variables)+size(model%surface_state_variables)))
call fabm_do_bottom(model,pelagic_rates(1:size(model%state_variables)), &
if (int2logical(do_bottom)) call fabm_do_bottom(model,pelagic_rates(1:size(model%state_variables)), &
pelagic_rates(size(model%state_variables)+size(model%surface_state_variables)+1:))
pelagic_rates(1:size(model%state_variables)) = pelagic_rates(1:size(model%state_variables))/column_depth
if (int2logical(do_surface) .or. int2logical(do_bottom)) then
if (.not.associated(column_depth)) call driver%fatal_error('get_rates', &
'Value for environmental dependency '//trim(environment_names(index_column_depth))// &
' must be provided if get_rates is called with the do_surface and/or do_bottom flags.')
pelagic_rates(1:size(model%state_variables)) = pelagic_rates(1:size(model%state_variables))/column_depth
end if
call fabm_do(model,pelagic_rates(1:size(model%state_variables)))
! Compute rate of change in conserved quantities
......
......@@ -24,7 +24,7 @@ except ImportError:
print 'Unable to load pyfabm. See https://github.com/fabm-model/code/wiki/python.'
sys.exit(1)
def evaluate(yaml_path, sources=(), location={}, assignments={}, verbose=True, ignore_missing=False):
def evaluate(yaml_path, sources=(), location={}, assignments={}, verbose=True, ignore_missing=False, surface=True, bottom=True):
# Create model object from YAML file.
model = pyfabm.Model(yaml_path)
......@@ -113,7 +113,7 @@ def evaluate(yaml_path, sources=(), location={}, assignments={}, verbose=True, i
print ' %s: %s %s' % (variable.name, variable.value, variable.units)
# Get model rates
rates = model.getRates()
rates = model.getRates(surface=surface, bottom=bottom)
assert len(rates) == len(model.state_variables), 'Length of array with rates does not match number of state variables'
if verbose:
......@@ -153,10 +153,12 @@ if __name__ == '__main__':
parser.add_argument('-l', '--location', nargs='*', help='NetCDF dimension to fix at particular index (specify: DIMENSION_NAME=INDEX)', default=[])
parser.add_argument('-v', '--values', nargs='*', help='Additional state variable/environmental dependency values (specify: VARIABLE_NAME=VALUE)', default=[])
parser.add_argument('--ignore_missing', action='store_true', help='Whether to ignore missing values for state variables and dependencies (the model will be evaluated with a default value of 0 for such missing variables)', default=False)
parser.add_argument('--no_surface', dest='surface', action='store_false', help='Whether to omit surface processes (do_surface calls)', default=True)
parser.add_argument('--no_bottom', dest='bottom', action='store_false', help='Whether to omit surface processes (do_bottom calls)', default=True)
parser.add_argument('--pause', action='store_true', help='Whether to pause before model evaluation to manually attach a debugger.', default=False)
args = parser.parse_args()
if args.pause:
raw_input('Attach the debugger (process id = %i) and then press Enter.' % os.getpid())
evaluate(args.model_path, args.sources, location=dict([dimension2index.split('=') for dimension2index in args.location]), assignments=dict([name2value.split('=') for name2value in args.values]), ignore_missing=args.ignore_missing)
evaluate(args.model_path, args.sources, location=dict([dimension2index.split('=') for dimension2index in args.location]), assignments=dict([name2value.split('=') for name2value in args.values]), ignore_missing=args.ignore_missing, surface=args.surface, bottom=args.bottom)
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