built-in models: add tracer and test case (#103)

* add builtin/tracer

* correct output visibility of result of single-term sums

* set units for result of aggregate sums and source sums

* add test case for builtin models

* use standard unit notation for sources

src/CMakeLists.txt

@@ -116,6 +116,7 @@ add_library(fabm_base OBJECT
             ${CMAKE_CURRENT_LIST_DIR}/fabm_types.F90
             ${CMAKE_CURRENT_LIST_DIR}/fabm_particle.F90
             ${CMAKE_CURRENT_LIST_DIR}/fabm_expressions.F90
+            ${CMAKE_CURRENT_LIST_DIR}/builtin/tracer.F90
             ${CMAKE_CURRENT_LIST_DIR}/builtin/constant.F90
             ${CMAKE_CURRENT_LIST_DIR}/builtin/reduction.F90
             ${CMAKE_CURRENT_LIST_DIR}/builtin/sum.F90

src/builtin/model_library.F90

@@ -10,6 +10,7 @@ module fabm_builtin_models
    use fabm_builtin_source
    use fabm_builtin_relaxation
    use fabm_builtin_depth_mapping
+   use fabm_builtin_tracer
 
    implicit none
 
@@ -62,6 +63,7 @@ subroutine create(self,name,model)
       class (type_base_model),pointer :: model
 
       select case (name)
+         case ('tracer');                   allocate(type_tracer::model)
          case ('bulk_constant');            allocate(type_interior_constant::model)
          case ('interior_constant');        allocate(type_interior_constant::model)
          case ('horizontal_constant');      allocate(type_horizontal_constant::model)
@@ -88,7 +90,7 @@ subroutine create(self,name,model)
          case ('bottom_temporal_mean');     allocate(type_bottom_temporal_mean::model)
          case ('surface_temporal_maximum'); allocate(type_surface_temporal_maximum::model)
          case ('depth_integrated_particle_override'); allocate(type_depth_integrated_particle_override::model)
-         case ('vertical_depth_range'); allocate(type_vertical_depth_range::model)
+         case ('vertical_depth_range');     allocate(type_vertical_depth_range::model)
          ! Add new examples models here
       end select
 

src/builtin/sum.F90

@@ -312,7 +312,7 @@ subroutine weighted_sum_initialize(self, configunit)
       integer,                   intent(in)            :: configunit
 
       type (type_component), pointer :: component
-      integer           :: i, n
+      integer           :: i, n, output
       character(len=10) :: temp
       class (type_weighted_sum_sms_distributor), pointer :: sms_distributor
       class (type_scaled_interior_variable), pointer :: scaled_variable
@@ -338,8 +338,10 @@ subroutine weighted_sum_initialize(self, configunit)
             output=self%result_output, source=source_constant, link=self%result_link)
          return
       elseif (n == 1) then
+         output = self%result_output
+         if (output /= output_none) output = ior(output, output_always_available)
          call self%add_interior_variable('result', self%units, 'result', link=self%result_link, &
-            act_as_state_variable=self%act_as_state_variable, presence=presence_external_required)
+            output=output, act_as_state_variable=self%act_as_state_variable, presence=presence_external_required)
          if (self%first%weight == 1.0_rk) then
             ! One component with scale factor 1 - directly link to the component's source variable.
             call request_coupling_to_component(self, self%result_link, self%first)
@@ -435,7 +437,7 @@ subroutine horizontal_weighted_sum_initialize(self, configunit)
       integer,                              intent(in)            :: configunit
 
       type (type_component), pointer :: component
-      integer           :: i, n
+      integer           :: i, n, output
       character(len=10) :: temp
       class (type_scaled_horizontal_variable), pointer :: scaled_variable
 
@@ -461,8 +463,10 @@ subroutine horizontal_weighted_sum_initialize(self, configunit)
          return
       elseif (n == 1) then
          ! One component only.
-         call self%add_horizontal_variable('result', self%units, 'result', output=self%result_output, link=self%result_link, domain=self%domain, &
-            act_as_state_variable=self%act_as_state_variable, presence=presence_external_required)
+         output = self%result_output
+         if (output /= output_none) output = ior(output, output_always_available)
+         call self%add_horizontal_variable('result', self%units, 'result', link=self%result_link, domain=self%domain, &
+            output=output, act_as_state_variable=self%act_as_state_variable, presence=presence_external_required)
          if (self%first%weight == 1.0_rk) then
             ! One component with scale factor 1 - directly link to the component's source variable.
             call request_coupling_to_component(self, self%result_link, self%first)

src/builtin/tracer.F90

@@ -0,0 +1,55 @@
+#include "fabm_driver.h"
+
+module fabm_builtin_tracer
+
+   ! This model describe a single passive tracer. Optionally, a vertical velocity
+   ! (sinking/floating) and light attenuation coefficient can be specified.
+   ! The unit is mol/m^3 by default, but may be explicitly configured as well.
+
+   use fabm_types
+
+   implicit none
+
+   private
+
+   type, extends(type_base_model), public :: type_tracer
+      ! Variables
+      type (type_state_variable_id) :: id_c
+   contains
+      procedure initialize
+   end type
+
+contains
+
+   subroutine initialize(self, configunit)
+      class (type_tracer), intent(inout), target :: self
+      integer,             intent(in)            :: configunit
+
+      character(len=64)               :: units
+      real(rk)                        :: vertical_velocity
+      real(rk)                        :: specific_light_attenuation
+      real(rk), parameter             :: days_per_second = 1.0_rk/86400.0_rk
+      logical                         :: conserved
+      character(len=attribute_length) :: standard_name
+
+      call self%register_implemented_routines()
+
+      ! Retrieve parameter values
+      call self%get_parameter(units, 'units', default='mol m-3')
+      call self%get_parameter(vertical_velocity, 'vertical_velocity', 'm d-1', 'vertical velocity (negative for settling, positive for rising)', default=0.0_rk, scale_factor=days_per_second)
+      call self%get_parameter(specific_light_attenuation, 'specific_light_attenuation', 'm-1 ('//trim(units)//')-1', 'specific light attenuation', default=0.0_rk)
+      call self%get_parameter(conserved, 'conserved', '', 'treat tracer as conserved quantity (activates budget tracking in host)', default=.false.)
+
+      ! Register state variables
+      call self%register_state_variable(self%id_c, 'c', units, 'concentration', 1.0_rk, &
+         minimum=0.0_rk, vertical_movement=vertical_velocity, specific_light_extinction=specific_light_attenuation)
+
+      if (conserved) then
+         ! User has activated conservation checking.
+         ! We create a "conserved quantity" for this tracer, which should induce the host model to compute and save the integral across the domain.
+         standard_name = get_safe_name(trim(self%get_path()) // '_total')
+         call self%add_to_aggregate_variable(type_universal_standard_variable(name=standard_name(2:), units=units, aggregate_variable=.true., conserved=.true.), self%id_c)
+      end if
+   end subroutine initialize
+
+end module

src/fabm_coupling.F90

@@ -473,6 +473,7 @@ subroutine create_sum(link, link_list, domain)
          end if
          sum%result_output = output_always_available
          sum%missing_value = 0.0_rk
+         sum%units = link%target%units
          component_link => link_list%first
          do while (associated(component_link))
             call sum%add_component(component_link)
@@ -649,6 +650,7 @@ recursive subroutine create_aggregate_models(self)
             horizontal_sum%domain = standard_variable2domain(standard_variable)
             sum => horizontal_sum
          end select
+         sum%units = aggregate_variable_access%link%target%units
          sum%act_as_state_variable = aggregate_variable_access%link%target%fake_state_variable
          sum%result_output = output_none
 

src/fabm_types.F90

@@ -1677,7 +1677,7 @@ subroutine register_source(self, link, sms_id, source)
       if (present(source)) source_ = source
       if (.not. self%implements(source_)) source_ = source_constant
       if (.not. associated(sms_id%link)) call self%add_interior_variable(trim(link%name)//'_sms', &
-         trim(link%target%units)//'/s', trim(link%target%long_name)//' sources-sinks', fill_value=0.0_rk, &
+         trim(link%target%units)//' s-1', trim(link%target%long_name)//' sources-sinks', fill_value=0.0_rk, &
          missing_value=0.0_rk, output=output_none, write_index=sms_id%sum_index, source=source_, link=sms_id%link)
       sms_id%link%target%write_operator = operator_add
       link2 => link%target%sms_list%append(sms_id%link%target, sms_id%link%target%name)
@@ -1697,7 +1697,7 @@ subroutine register_surface_flux(self, link, surface_flux_id, source)
       if (present(source)) source_ = source
       if (.not. self%implements(source_)) source_ = source_constant
       if (.not. associated(surface_flux_id%link)) call self%add_horizontal_variable(trim(link%name) // '_sfl', &
-         trim(link%target%units) // '*m/s', trim(link%target%long_name) // ' surface flux', fill_value=0.0_rk, &
+         trim(link%target%units) // ' m s-1', trim(link%target%long_name) // ' surface flux', fill_value=0.0_rk, &
          missing_value=0.0_rk, output=output_none, write_index=surface_flux_id%horizontal_sum_index, &
          domain=domain_surface, source=source_, link=surface_flux_id%link)
       surface_flux_id%link%target%write_operator = operator_add
@@ -1718,7 +1718,7 @@ subroutine register_bottom_flux(self, link, bottom_flux_id, source)
       if (present(source)) source_ = source
       if (.not. self%implements(source_)) source_ = source_constant
       if (.not. associated(bottom_flux_id%link)) call self%add_horizontal_variable(trim(link%name) // '_bfl', &
-         trim(link%target%units) // '*m/s', trim(link%target%long_name) // ' bottom flux', fill_value=0.0_rk, &
+         trim(link%target%units) // ' m s-1', trim(link%target%long_name) // ' bottom flux', fill_value=0.0_rk, &
          missing_value=0.0_rk, output=output_none, write_index=bottom_flux_id%horizontal_sum_index, &
          domain=domain_bottom, source=source_, link=bottom_flux_id%link)
       bottom_flux_id%link%target%write_operator = operator_add
@@ -1738,7 +1738,7 @@ subroutine register_movement(self, link, movement_id, vertical_movement)
       vertical_movement_ = 0.0_rk
       if (present(vertical_movement)) vertical_movement_ = vertical_movement
       if (.not. associated(movement_id%link)) call self%add_interior_variable(trim(link%name) // '_w', &
-         'm/s', trim(link%target%long_name) // ' vertical velocity', fill_value=vertical_movement_, missing_value=0.0_rk, &
+         'm s-1', trim(link%target%long_name) // ' vertical velocity', fill_value=vertical_movement_, missing_value=0.0_rk, &
          output=output_none, write_index=movement_id%sum_index, link=movement_id%link, source=source_constant)
       if (self%implements(source_get_vertical_movement)) then
          movement_id%link%target%source = source_get_vertical_movement
@@ -1760,7 +1760,7 @@ subroutine register_surface_source(self, link, sms_id, source)
       if (present(source)) source_ = source
       if (.not. self%implements(source_)) source_ = source_constant
       if (.not. associated(sms_id%link)) call self%add_horizontal_variable(trim(link%name) // '_sms', &
-         trim(link%target%units) // '/s', trim(link%target%long_name) // ' sources-sinks', fill_value=0.0_rk, &
+         trim(link%target%units) // ' s-1', trim(link%target%long_name) // ' sources-sinks', fill_value=0.0_rk, &
          missing_value=0.0_rk, output=output_none, write_index=sms_id%horizontal_sum_index, link=sms_id%link, &
          domain=domain_surface, source=source_)
       sms_id%link%target%write_operator = operator_add
@@ -1781,7 +1781,7 @@ subroutine register_bottom_source(self, link, sms_id, source)
       if (present(source)) source_ = source
       if (.not. self%implements(source_)) source_ = source_constant
       if (.not. associated(sms_id%link)) call self%add_horizontal_variable(trim(link%name) // '_sms', &
-         trim(link%target%units) // '/s', trim(link%target%long_name) // ' sources-sinks', fill_value=0.0_rk, &
+         trim(link%target%units) // ' s-1', trim(link%target%long_name) // ' sources-sinks', fill_value=0.0_rk, &
          missing_value=0.0_rk, output=output_none, write_index=sms_id%horizontal_sum_index, link=sms_id%link, &
          domain=domain_bottom, source=source_)
       sms_id%link%target%write_operator = operator_add

testcases/fabm-builtin.yaml

@@ -0,0 +1,126 @@
+instances:
+  tracer:
+    model: tracer
+  interior_constant:
+    model: interior_constant
+    parameters:
+      value: 1.0
+  horizontal_constant:
+    model: horizontal_constant
+    parameters:
+      value: 1.0
+  surface_constant:
+    model: horizontal_constant
+    parameters:
+      value: 1.0
+  bottom_constant:
+    model: horizontal_constant
+    parameters:
+      value: 1.0
+  global_constant:
+    model: global_constant
+    parameters:
+      value: 1.0
+  interior_sum_empty:
+    model: weighted_sum
+    parameters:
+      n: 0
+  interior_sum_1term:
+    model: weighted_sum
+    parameters:
+      n: 1
+    coupling:
+      term1: interior_constant/data
+  interior_sum_nterms:
+    model: weighted_sum
+    parameters:
+      n: 2
+    coupling:
+      term1: tracer/c
+      term2: interior_constant/data
+  interior_sum_1term_scaled:
+    model: weighted_sum
+    parameters:
+      n: 1
+      weight1: 2.0
+    coupling:
+      term1: interior_constant/data
+  interior_sum_nterm_scaled:
+    model: weighted_sum
+    parameters:
+      n: 2
+      weight1: 1.0
+      weight2: 2.0
+    coupling:
+      term1: tracer/c
+      term2: interior_constant/data      
+  horizontal_sum_empty:
+    model: horizontal_weighted_sum
+    parameters:
+      n: 0
+  horizontal_sum_1term:
+    model: horizontal_weighted_sum
+    parameters:
+      n: 1
+    coupling:
+      term1: horizontal_constant/data
+  horizontal_sum_nterms:
+    model: horizontal_weighted_sum
+    parameters:
+      n: 3
+    coupling:
+      term1: horizontal_constant/data
+      term2: surface_constant/data
+      term3: bottom_constant/data      
+  horizontal_sum_1term_scaled:
+    model: horizontal_weighted_sum
+    parameters:
+      n: 1
+      weight1: 2.0
+    coupling:
+      term1: horizontal_constant/data
+  horizontal_sum_nterms_scaled:
+    model: horizontal_weighted_sum
+    parameters:
+      n: 3
+      weight1: 1.0
+      weight2: 2.0
+      weight3: 3.0            
+    coupling:
+      term1: horizontal_constant/data
+      term2: surface_constant/data
+      term3: bottom_constant/data
+  bottom_layer:
+    model: bottom_layer
+    coupling:
+      source: tracer/c
+  surface_layer:
+    model: surface_layer
+    coupling:
+      source: tracer/c
+  vertical_integral:
+    model: vertical_integral
+    coupling:
+      source: tracer/c
+  vertical_average:
+    model: vertical_integral
+    parameters:
+      average: true
+    coupling:
+      source: tracer/c
+  bounded_vertical_integral:
+    model: bounded_vertical_integral
+    parameters:
+      minimum_depth: 10.0
+      maximum_depth: 20.0
+    coupling:
+      source: tracer/c
+  bounded_vertical_average:
+    model: bounded_vertical_integral
+    parameters:
+      average: true
+      minimum_depth: 10.0
+      maximum_depth: 20.0
+    coupling:
+      source: tracer/c
+