pyfabm: simplified cell_thickness treatment

src/c/fabm_c.F90

@@ -37,12 +37,10 @@ module fabm_c
    end type
 
    type type_model_wrapper
-      class (type_fabm_model), pointer           :: p => null()
-      type (type_variable_list)                  :: environment
-      type (type_internal_variable), pointer     :: cell_thickness_variable
-      type (type_link_list)                      :: coupling_link_list
-      real(c_double) _DIMENSION_GLOBAL_, pointer :: cell_thickness => null()
-      type (type_property_dictionary)            :: forced_parameters, forced_couplings
+      class (type_fabm_model), pointer :: p => null()
+      type (type_variable_list)        :: environment
+      type (type_link_list)            :: coupling_link_list
+      type (type_property_dictionary)  :: forced_parameters, forced_couplings
    end type
 
 contains
@@ -96,7 +94,7 @@ function create_model(path _POSTARG_LOCATION_) result(ptr) bind(c)
       call model%p%set_domain(_PREARG_LOCATION_ 1._rk)
 
       ! Retrieve arrays to hold values for environmental variables and corresponding metadata.
-      call get_environment_metadata(model%p, model%environment, model%cell_thickness_variable)
+      call get_environment_metadata(model%p, model%environment)
 
       call get_couplings(model%p, model%coupling_link_list)
 
@@ -163,8 +161,7 @@ subroutine reinitialize(model)
       call model%p%set_domain(_PREARG_LOCATION_ 1._rk)
 
       ! Retrieve arrays to hold values for environmental variables and corresponding metadata.
-      call get_environment_metadata(model%p, model%environment, model%cell_thickness_variable)
-      model%cell_thickness => null()
+      call get_environment_metadata(model%p, model%environment)
 
       call get_couplings(model%p, model%coupling_link_list)
    end subroutine reinitialize
@@ -431,16 +428,12 @@ subroutine link_dependency_data(pmodel, pvariable, dat) bind(c)
       case (domain_interior)
          interior_data => c_f_pointer_interior(model, dat)
          call model%p%link_interior_data(variable, interior_data)
-         if (associated(variable, model%cell_thickness_variable)) model%cell_thickness => interior_data
       case (domain_scalar)
          call c_f_pointer(c_loc(dat), scalar_data)
          call model%p%link_scalar(variable, scalar_data)
       case default
          horizontal_data => c_f_pointer_horizontal(model, dat)
          call model%p%link_horizontal_data(variable, horizontal_data)
-#ifndef _FABM_DEPTH_DIMENSION_INDEX_
-         if (associated(variable, model%cell_thickness_variable)) model%cell_thickness => horizontal_data
-#endif
       end select
    end subroutine link_dependency_data
 
@@ -489,19 +482,20 @@ subroutine link_bottom_state_data(pmodel, index, dat) bind(c)
       call model%p%link_bottom_state_data(index, dat_)
    end subroutine link_bottom_state_data
 
-   subroutine get_sources(pmodel, t, sources_interior, sources_surface, sources_bottom, do_surface, do_bottom) bind(c)
+   subroutine get_sources(pmodel, t, sources_interior, sources_surface, sources_bottom, do_surface, do_bottom, cell_thickness) bind(c)
       !DIR$ ATTRIBUTES DLLEXPORT :: get_sources
       type (c_ptr),   value,  intent(in) :: pmodel
       real(c_double), value,  intent(in) :: t
       real(c_double), target, intent(in) :: sources_interior(*), sources_surface(*), sources_bottom(*)
       integer(c_int), value,  intent(in) :: do_surface, do_bottom
+      real(c_double), target, intent(in) :: cell_thickness(*)
 
       logical :: surface, bottom
       type (type_model_wrapper), pointer :: model
       real(c_double) _DIMENSION_GLOBAL_PLUS_1_, pointer :: sources_interior_
       real(c_double) _DIMENSION_GLOBAL_HORIZONTAL_PLUS_1_, pointer :: sources_surface_, sources_bottom_
       real(c_double) _DIMENSION_HORIZONTAL_SLICE_PLUS_1_, allocatable :: fluxes
-      real(c_double) :: cell_thickness
+      real(c_double) _ATTRIBUTES_GLOBAL_, pointer :: cell_thickness_
       _DECLARE_LOCATION_
 
 #  if _FABM_DIMENSION_COUNT_ > 0
@@ -529,17 +523,13 @@ subroutine get_sources(pmodel, t, sources_interior, sources_surface, sources_bot
 
       surface = int2logical(do_surface)
       bottom = int2logical(do_bottom)
-      if ((surface .or. bottom) .and. size(model%p%interior_state_variables) > 0 .and. .not. associated(model%cell_thickness)) then
-          call driver%fatal_error('get_sources', &
-            'Value for cell_thickness must be provided if get_sources is called with the do_surface and/or do_bottom flags.')
-          return
-      end if
+      if ((surface .or. bottom) .and. size(model%p%interior_state_variables) > 0) cell_thickness_ => c_f_pointer_interior(model, cell_thickness)
 
       call c_f_pointer(c_loc(sources_interior), sources_interior_, (/_PREARG_LOCATION_ size(model%p%interior_state_variables)/))
       call c_f_pointer(c_loc(sources_surface), sources_surface_, (/_PREARG_HORIZONTAL_LOCATION_ size(model%p%surface_state_variables)/))
       call c_f_pointer(c_loc(sources_bottom), sources_bottom_, (/_PREARG_HORIZONTAL_LOCATION_ size(model%p%bottom_state_variables)/))
 #ifdef _HORIZONTAL_IS_VECTORIZED_
-      allocate(fluxes(_ITERATOR_,size(model%p%interior_state_variables)))
+      allocate(fluxes(_ITERATOR_, size(model%p%interior_state_variables)))
 #else
       allocate(fluxes(size(model%p%interior_state_variables)))
 #endif
@@ -574,11 +564,11 @@ subroutine get_sources(pmodel, t, sources_interior, sources_surface, sources_bot
 #ifdef _HORIZONTAL_IS_VECTORIZED_
                _DO_CONCURRENT_(_ITERATOR_,_START_,_STOP_)
                   sources_interior_(_PREARG_LOCATION_ :) = sources_interior_(_PREARG_LOCATION_ :) &
-                     + fluxes(_ITERATOR_,:) / model%cell_thickness _INDEX_LOCATION_
+                     + fluxes(_ITERATOR_,:) / cell_thickness_ _INDEX_LOCATION_
                end do
 #else
                sources_interior_(_PREARG_LOCATION_ :) = sources_interior_(_PREARG_LOCATION_ :) &
-                  + fluxes(:) / model%cell_thickness _INDEX_LOCATION_
+                  + fluxes(:) / cell_thickness_ _INDEX_LOCATION_
 #endif
             end if
          _END_OUTER_HORIZONTAL_LOOP_
@@ -607,14 +597,14 @@ subroutine get_sources(pmodel, t, sources_interior, sources_surface, sources_bot
                   _VERTICAL_ITERATOR_ = model%p%domain%bottom_indices _INDEX_HORIZONTAL_LOCATION_
 #endif
                   sources_interior_(_PREARG_LOCATION_ :) = sources_interior_(_PREARG_LOCATION_ :) &
-                     + fluxes(_ITERATOR_,:) / model%cell_thickness _INDEX_LOCATION_
+                     + fluxes(_ITERATOR_,:) / cell_thickness_ _INDEX_LOCATION_
                end do
 #else
 #if _FABM_BOTTOM_INDEX_==-1
                _VERTICAL_ITERATOR_ = model%p%domain%bottom_indices _INDEX_HORIZONTAL_LOCATION_
 #endif
                sources_interior_(_PREARG_LOCATION_ :) = sources_interior_(_PREARG_LOCATION_ :) &
-                  + fluxes(:) / model%cell_thickness _INDEX_LOCATION_
+                  + fluxes(:) / cell_thickness_ _INDEX_LOCATION_
 #endif
             end if
          _END_OUTER_HORIZONTAL_LOOP_

src/c/helper.F90

@@ -13,15 +13,12 @@ module fabm_python_helper
 
 contains
 
-   subroutine get_environment_metadata(model, list, cell_thickness)
+   subroutine get_environment_metadata(model, list)
      type (type_fabm_model),    intent(inout) :: model
      type (type_variable_list), intent(out)   :: list
-     type (type_internal_variable), pointer   :: cell_thickness
 
      type (type_link), pointer :: link
 
-     cell_thickness => null()
-
       ! Get number of environmental dependencies (light, temperature, etc.)
       link => model%links_postcoupling%first
       do while (associated(link))
@@ -30,16 +27,8 @@ subroutine get_environment_metadata(model, list, cell_thickness)
             select case (link%target%domain)
             case (domain_interior)
                if (.not. associated(model%catalog%interior(link%target%catalog_index)%p)) call list%append(link%target)
-               if (.not. associated(cell_thickness) .and. associated(link%target%standard_variables%first)) then
-                  if (link%target%standard_variables%contains(standard_variables%cell_thickness)) cell_thickness => link%target
-               end if
             case (domain_bottom, domain_surface, domain_horizontal)
                if (.not. associated(model%catalog%horizontal(link%target%catalog_index)%p)) call list%append(link%target)
-#ifndef _FABM_DEPTH_DIMENSION_INDEX_
-               if (.not. associated(cell_thickness) .and. associated(link%target%standard_variables%first)) then
-                  if (link%target%standard_variables%contains(standard_variables%bottom_depth)) cell_thickness => link%target
-               end if
-#endif
             case (domain_scalar)
                if (.not. associated(model%catalog%scalar(link%target%catalog_index)%p)) call list%append(link%target)
             end select

src/c/integrate.F90

@@ -1,4 +1,5 @@
 #include "fabm_driver.h"
+#include "fabm_private.h"
 
 module fabm_c_integrate
 #if _FABM_DIMENSION_COUNT_ == 0
@@ -13,13 +14,14 @@ module fabm_c_integrate
 
 contains
 
-   subroutine integrate(pmodel, nt, ny, t_, y_ini_, y_, dt, do_surface, do_bottom) bind(c)
+   subroutine integrate(pmodel, nt, ny, t_, y_ini_, y_, dt, do_surface, do_bottom, cell_thickness) bind(c)
       !DIR$ ATTRIBUTES DLLEXPORT :: integrate
       type (c_ptr),  value, intent(in) :: pmodel
       integer(c_int),value, intent(in) :: nt, ny
       real(c_double),target,intent(in) :: t_(*), y_ini_(*), y_(*)
       real(c_double),value, intent(in) :: dt
       integer(c_int),value, intent(in) :: do_surface, do_bottom
+      real(c_double), target, intent(in) :: cell_thickness(*)
 
       type (type_model_wrapper), pointer :: model
       real(c_double), pointer :: t(:), y_ini(:), y(:,:)
@@ -28,6 +30,7 @@ subroutine integrate(pmodel, nt, ny, t_, y_ini_, y_, dt, do_surface, do_bottom)
       real(rke), target      :: y_cur(ny)
       real(rke)              :: dy(ny)
       logical                :: surface, bottom
+      real(c_double) _ATTRIBUTES_GLOBAL_, pointer :: cell_thickness_
 
       call c_f_pointer(pmodel, model)
       if (model%p%status < status_start_done) then
@@ -46,11 +49,7 @@ subroutine integrate(pmodel, nt, ny, t_, y_ini_, y_, dt, do_surface, do_bottom)
 
       surface = int2logical(do_surface)
       bottom = int2logical(do_bottom)
-      if ((surface .or. bottom) .and. .not. associated(model%cell_thickness)) then
-          call driver%fatal_error('integrate', &
-            'Value for cell_thickness must be provided if integrate is called with the do_surface and/or do_bottom flags.')
-          return
-      end if
+      cell_thickness_ => c_f_pointer_interior(model, cell_thickness)
       call model%p%link_all_interior_state_data(y_cur(1:size(model%p%interior_state_variables)))
       call model%p%link_all_surface_state_data(y_cur(size(model%p%interior_state_variables) + 1: &
          size(model%p%interior_state_variables) + size(model%p%surface_state_variables)))
@@ -74,7 +73,7 @@ subroutine integrate(pmodel, nt, ny, t_, y_ini_, y_, dt, do_surface, do_bottom)
           if (bottom) call model%p%get_bottom_sources(dy(1:size(model%p%interior_state_variables)), &
              dy(size(model%p%interior_state_variables) + size(model%p%surface_state_variables) + 1:))
           if (surface .or. bottom) dy(1:size(model%p%interior_state_variables)) = dy(1:size(model%p%interior_state_variables)) &
-             / model%cell_thickness
+             / cell_thickness_
           call model%p%get_interior_sources(dy(1:size(model%p%interior_state_variables)))
           y_cur = y_cur + dt * dy * 86400
           t_cur = t_cur + dt

src/drivers/python/pyfabm/__init__.py

@@ -135,11 +135,11 @@ def find_library(basedir):
 fabm.start.restype = None
 
 # Routine for retrieving source-sink terms for the interior domain.
-fabm.get_sources.argtypes = [ctypes.c_void_p, ctypes.c_double, numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), ctypes.c_int, ctypes.c_int]
+fabm.get_sources.argtypes = [ctypes.c_void_p, ctypes.c_double, numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)]
 fabm.get_sources.restype = None
 fabm.check_state.argtypes = [ctypes.c_void_p, ctypes.c_int]
 fabm.check_state.restype = ctypes.c_int
-fabm.integrate.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int, numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags=CONTIGUOUS), ctypes.c_double, ctypes.c_int, ctypes.c_int]
+fabm.integrate.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int, numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1, flags=CONTIGUOUS), numpy.ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags=CONTIGUOUS), ctypes.c_double, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)]
 fabm.integrate.restype = None
 
 # Routine for getting git repository version information.
@@ -186,6 +186,9 @@ def reploldminus(m):
     #unit = oldsupminus.sub(reploldminus,unit)
     return unit
 
+class FABMException(Exception):
+    pass
+
 def hasError():
    return fabm.get_error_state() != 0
 
@@ -399,13 +402,22 @@ def __init__(self, pmodel, name):
         self.user_created = iuser.value!=0
 
 class Model(object):
-    def __init__(self,path='fabm.yaml'):
+    def __init__(self, path='fabm.yaml'):
         fabm.reset_error_state()
         self.lookup_tables = {}
+        self._cell_thickness = None
         self.pmodel = fabm.create_model(path.encode('ascii'))
-        assert not hasError(), 'An error occurred while parsing %s:\n%s' % (path, getError())
+        if hasError():
+            raise FABMException('An error occurred while parsing %s:\n%s' % (path, getError()))
         self.updateConfiguration()
 
+    def setCellThickness(self, value):
+        if self._cell_thickness is None:
+            self._cell_thickness = ctypes.c_double()
+        self._cell_thickness.value = value
+
+    cell_thickness = property(fset=setCellThickness)
+
     def getSubModel(self,name):
         return SubModel(self.pmodel, name)
 
@@ -497,7 +509,8 @@ def updateConfiguration(self,settings=None):
         self.state_variables = self.interior_state_variables + self.surface_state_variables + self.bottom_state_variables
         self.diagnostic_variables = self.interior_diagnostic_variables + self.horizontal_diagnostic_variables
 
-        if settings is not None: self.restoreSettings(settings)
+        if settings is not None:
+            self.restoreSettings(settings)
 
         # For backward compatibility
         self.bulk_state_variables = self.interior_state_variables
@@ -515,11 +528,17 @@ def getRates(self, t=None, surface=True, bottom=True):
         sources_interior = sources[:len(self.interior_state_variables)]
         sources_surface = sources[len(self.interior_state_variables):len(self.interior_state_variables)+len(self.surface_state_variables)]
         sources_bottom = sources[len(self.interior_state_variables)+len(self.surface_state_variables):]
-        fabm.get_sources(self.pmodel, t, sources_interior, sources_surface, sources_bottom, surface, bottom)
+        assert not ((surface or bottom) and self._cell_thickness is None), 'You must assign model.cell_thickness to use getRates'
+        fabm.get_sources(self.pmodel, t, sources_interior, sources_surface, sources_bottom, surface, bottom, ctypes.byref(self._cell_thickness))
+        if hasError():
+            raise FABMException(getError())
         return sources
 
     def checkState(self, repair=False):
-        return fabm.check_state(self.pmodel, repair) != 0
+        valid = fabm.check_state(self.pmodel, repair) != 0
+        if hasError():
+            raise FABMException(getError())
+        return valid
 
     def getJacobian(self,pert=None):
         # Define perturbation per state variable.
@@ -635,11 +654,14 @@ def printArray(classname, array):
 
 class Simulator(object):
     def __init__(self, model):
+        assert model._cell_thickness is not None, 'You must assign model.cell_thickness to use Simulator'
         self.model = model
 
     def integrate(self, y0, t, dt, surface=True, bottom=True):
         y = numpy.empty((t.size, self.model.state.size))
-        fabm.integrate(self.model.pmodel, t.size, self.model.state.size, t, y0, y, dt, surface, bottom)
+        fabm.integrate(self.model.pmodel, t.size, self.model.state.size, t, y0, y, dt, surface, bottom, ctypes.byref(self.model._cell_thickness))
+        if hasError():
+            raise FABMException(getError())
         return y
 
 def unload():

util/developers/run_all_testcases.py

@@ -192,6 +192,7 @@ def test_pyfabm(args, testcases):
         print('  %s... ' % case, end='')
         sys.stdout.flush()
         m = pyfabm.Model(path)
+        m.cell_thickness = 1.
         for d in m.dependencies:
             dependency_names.add(d.name)
             d.value = 1.