diff --git a/CMakeLists.txt b/CMakeLists.txt index f618578a..5afa010b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,7 @@ cmake_minimum_required(VERSION 3.13) project(fabm - VERSION 2.1.6 + VERSION 3.0.0 DESCRIPTION "Framework for Aquatic Biogeochemical Models" HOMEPAGE_URL https://fabm.net LANGUAGES Fortran diff --git a/include/fabm.h b/include/fabm.h index dfdc498a..45c70e12 100644 --- a/include/fabm.h +++ b/include/fabm.h @@ -107,9 +107,9 @@ #else ! Interior procedures operate in 0D ! Interior slices may be 0D scalars or 1D arrays [the latter if the model has a vertical dimension] -# define _LOOP_BEGIN_EX_(cache) -# define _CONCURRENT_LOOP_BEGIN_EX_(cache) -# define _LOOP_END_ +# define _LOOP_BEGIN_EX_(cache) if (cache%n /= 0) then +# define _CONCURRENT_LOOP_BEGIN_EX_(cache) if (cache%n /= 0) then +# define _LOOP_END_ end if # ifdef _INTERIOR_IS_VECTORIZED_ ! Interior slices are 1D arrays - we will operate on their first element (_I_=1) # define _DECLARE_INTERIOR_INDICES_ integer,parameter :: _I_=1 @@ -134,9 +134,9 @@ #else ! Horizontal procedures operate in 0D ! Horizontal slices MUST be scalars; interior slices can be scalars or 1D arrays -# define _HORIZONTAL_LOOP_BEGIN_EX_(cache) -# define _CONCURRENT_HORIZONTAL_LOOP_BEGIN_EX_(cache) -# define _HORIZONTAL_LOOP_END_ +# define _HORIZONTAL_LOOP_BEGIN_EX_(cache) if (cache%n /= 0) then +# define _CONCURRENT_HORIZONTAL_LOOP_BEGIN_EX_(cache) if (cache%n /= 0) then +# define _HORIZONTAL_LOOP_END_ end if # ifdef _INTERIOR_IS_VECTORIZED_ ! Interior slices are 1D arrays - we will operate on their first element (_I_=1) # define _DECLARE_HORIZONTAL_INDICES_ integer,parameter :: _I_=1 @@ -186,11 +186,11 @@ ! Vertical procedures operate in 0D ! Interior slices may scalars or 1D arrays [the latter if the model is vectorized over a horizontal dimension] ! Applies to all models without depth dimension; for instance, 0D box or model with i,j or i -# define _CONCURRENT_VERTICAL_LOOP_BEGIN_EX_(cache) -# define _VERTICAL_LOOP_END_ +# define _CONCURRENT_VERTICAL_LOOP_BEGIN_EX_(cache) if (cache%n /= 0) then +# define _VERTICAL_LOOP_END_ end if # define _VERTICAL_LOOP_EXIT_ -# define _DOWNWARD_LOOP_BEGIN_ -# define _UPWARD_LOOP_BEGIN_ +# define _DOWNWARD_LOOP_BEGIN_ if (cache%n /= 0) then +# define _UPWARD_LOOP_BEGIN_ if (cache%n /= 0) then # define _MOVE_TO_SURFACE_ # define _MOVE_TO_BOTTOM_ # ifdef _INTERIOR_IS_VECTORIZED_ diff --git a/include/fabm_private.h b/include/fabm_private.h index e69454b2..604b8201 100644 --- a/include/fabm_private.h +++ b/include/fabm_private.h @@ -269,12 +269,18 @@ # define _IS_UNMASKED_(maskvalue) .true. #endif -#if defined(_FABM_CONTIGUOUS_)&&!defined(_FABM_NO_CONTIGUOUS_) +#ifndef _FABM_NO_CONTIGUOUS_ # define _CONTIGUOUS_ ,contiguous #else # define _CONTIGUOUS_ #endif +#ifdef _FABM_CONTIGUOUS_ +# define _HOST_CONTIGUOUS_ _CONTIGUOUS_ +#else +# define _HOST_CONTIGUOUS_ +#endif + #ifndef _BEGIN_OUTER_HORIZONTAL_LOOP_ # define _BEGIN_OUTER_HORIZONTAL_LOOP_ # define _END_OUTER_HORIZONTAL_LOOP_ @@ -326,7 +332,8 @@ #if _HORIZONTAL_DIMENSION_COUNT_>0 # define _INDEX_HORIZONTAL_LOCATION_ (_HORIZONTAL_LOCATION_) # define _DIMENSION_GLOBAL_HORIZONTAL_ ,dimension(_HORIZONTAL_LOCATION_DIMENSIONS_) -# define _ATTRIBUTES_GLOBAL_HORIZONTAL_ _DIMENSION_GLOBAL_HORIZONTAL_ _CONTIGUOUS_ +# define _ATTRIBUTES_GLOBAL_HORIZONTAL_ _DIMENSION_GLOBAL_HORIZONTAL_ _HOST_CONTIGUOUS_ +# define _ATTRIBUTES_GLOBAL_HORIZONTAL_CONTIGUOUS_ _DIMENSION_GLOBAL_HORIZONTAL_ _CONTIGUOUS_ # define _ARG_HORIZONTAL_LOCATION_ _HORIZONTAL_LOCATION_ # define _POSTARG_HORIZONTAL_LOCATION_ ,_ARG_HORIZONTAL_LOCATION_ # define _POSTARG_HORIZONTAL_LOCATION_RANGE_ ,_HORIZONTAL_LOCATION_RANGE_ @@ -339,6 +346,7 @@ # define _INDEX_HORIZONTAL_LOCATION_ # define _DIMENSION_GLOBAL_HORIZONTAL_ # define _ATTRIBUTES_GLOBAL_HORIZONTAL_ +# define _ATTRIBUTES_GLOBAL_HORIZONTAL_CONTIGUOUS_ # define _ARG_HORIZONTAL_LOCATION_ # define _POSTARG_HORIZONTAL_LOCATION_ # define _POSTARG_HORIZONTAL_LOCATION_RANGE_ @@ -356,7 +364,8 @@ #if _FABM_DIMENSION_COUNT_>0 # define _INDEX_LOCATION_ (_LOCATION_) # define _DIMENSION_GLOBAL_ ,dimension(_LOCATION_DIMENSIONS_) -# define _ATTRIBUTES_GLOBAL_ _DIMENSION_GLOBAL_ _CONTIGUOUS_ +# define _ATTRIBUTES_GLOBAL_ _DIMENSION_GLOBAL_ _HOST_CONTIGUOUS_ +# define _ATTRIBUTES_GLOBAL_CONTIGUOUS_ _DIMENSION_GLOBAL_ _CONTIGUOUS_ # define _POSTARG_LOCATION_ ,_LOCATION_ # define _POSTARG_LOCATION_RANGE_ ,_LOCATION_RANGE_ # define _PREARG_LOCATION_ _LOCATION_, @@ -368,6 +377,7 @@ # define _INDEX_LOCATION_ # define _DIMENSION_GLOBAL_ # define _ATTRIBUTES_GLOBAL_ +# define _ATTRIBUTES_GLOBAL_CONTIGUOUS_ # define _POSTARG_LOCATION_ # define _POSTARG_LOCATION_RANGE_ # define _PREARG_LOCATION_ diff --git a/include/fabm_version.h b/include/fabm_version.h index 43e7dedb..cfa4c097 100644 --- a/include/fabm_version.h +++ b/include/fabm_version.h @@ -1 +1 @@ -#define _FABM_API_VERSION_ 2 +#define _FABM_API_VERSION_ 3 diff --git a/pyproject.toml b/pyproject.toml index cdf3caa3..37cf1fb7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "pyfabm" -version = "2.1.6" +version = "3.0.0" authors = [ {name = "Jorn Bruggeman", email = "jorn@bolding-bruggeman.com"}, {name = "Karsten Bolding", email = "karsten@bolding-bruggeman.com"} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 54c8a092..27f2b804 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -115,6 +115,7 @@ add_library(fabm_base OBJECT ${CMAKE_CURRENT_LIST_DIR}/fabm_standard_variables.F90 ${CMAKE_CURRENT_LIST_DIR}/fabm_properties.F90 ${CMAKE_CURRENT_LIST_DIR}/fabm_types.F90 + ${CMAKE_CURRENT_LIST_DIR}/fabm_global_types.F90 ${CMAKE_CURRENT_LIST_DIR}/fabm_particle.F90 ${CMAKE_CURRENT_LIST_DIR}/fabm_expressions.F90 ${CMAKE_CURRENT_LIST_DIR}/builtin/tracer.F90 diff --git a/src/builtin/constant.F90 b/src/builtin/constant.F90 index e6e77de2..d9b117e3 100644 --- a/src/builtin/constant.F90 +++ b/src/builtin/constant.F90 @@ -45,7 +45,6 @@ subroutine interior_constant_initialize(self, configunit) real(rk) :: value type (type_interior_standard_variable) :: standard_variable - call self%register_implemented_routines() call self%get_parameter(standard_variable%name, 'standard_name', '', 'standard name', default='') call self%get_parameter(value, 'value', '', 'value') if (standard_variable%name /= '') then @@ -66,7 +65,6 @@ subroutine horizontal_constant_initialize(self, configunit) real(rk) :: value type (type_horizontal_standard_variable) :: standard_variable - call self%register_implemented_routines() call self%get_parameter(standard_variable%name, 'standard_name', '', 'standard name', default='') call self%get_parameter(value, 'value', '', 'value') if (standard_variable%name /= '') then @@ -87,7 +85,6 @@ subroutine surface_constant_initialize(self, configunit) real(rk) :: value type (type_surface_standard_variable) :: standard_variable - call self%register_implemented_routines() call self%get_parameter(standard_variable%name, 'standard_name', '', 'standard name', default='') call self%get_parameter(value, 'value', '', 'value') if (standard_variable%name /= '') then @@ -108,7 +105,6 @@ subroutine bottom_constant_initialize(self, configunit) real(rk) :: value type (type_bottom_standard_variable) :: standard_variable - call self%register_implemented_routines() call self%get_parameter(standard_variable%name, 'standard_name', '', 'standard name', default='') call self%get_parameter(value, 'value', '', 'value') if (standard_variable%name /= '') then @@ -129,7 +125,6 @@ subroutine global_constant_initialize(self, configunit) real(rk) :: value type (type_global_standard_variable) :: standard_variable - call self%register_implemented_routines() call self%get_parameter(standard_variable%name, 'standard_name', '', 'standard name', default='') call self%get_parameter(value, 'value', '', 'value') if (standard_variable%name /= '') then diff --git a/src/builtin/depth_integral.F90 b/src/builtin/depth_integral.F90 index 7c01365a..61748eb7 100644 --- a/src/builtin/depth_integral.F90 +++ b/src/builtin/depth_integral.F90 @@ -34,7 +34,6 @@ subroutine depth_integral_initialize(self, configunit) class (type_depth_integral), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_column/)) call self%get_parameter(self%average, 'average', '', 'compute average by dividing integral by total height', default=self%average) call self%register_dependency(self%id_input, 'source', '', 'source') call self%register_dependency(self%id_thickness, standard_variables%cell_thickness) diff --git a/src/builtin/depth_mapping.F90 b/src/builtin/depth_mapping.F90 index f9fe3a1d..dc3d7fc7 100644 --- a/src/builtin/depth_mapping.F90 +++ b/src/builtin/depth_mapping.F90 @@ -153,7 +153,6 @@ subroutine vertical_depth_range_initialize(self, configunit) class (type_vertical_depth_range), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do/)) call self%get_parameter(self%minimum_depth, 'minimum_depth', 'm', 'minimum depth (distance from surface)', & default=self%minimum_depth) call self%get_parameter(self%maximum_depth, 'maximum_depth', 'm', 'maximum depth (distance from surface)', & @@ -258,8 +257,6 @@ subroutine weighted_depth_integral_initialize(self, configunit) class (type_absolute_rate_distributor), pointer :: rate_distributor - call self%register_implemented_routines((/source_do_column/)) - call self%get_parameter(self%average, 'average', '', 'average', default=self%average) call self%get_parameter(self%act_as_state_variable, 'act_as_state_variable', '', 'act as state variable', & default=self%act_as_state_variable) @@ -324,8 +321,6 @@ subroutine projector_initialize(self, configunit) class (type_projected_rate_distributor), pointer :: rate_distributor - call self%register_implemented_routines((/source_do/)) - call self%get_parameter(self%act_as_state_variable, 'act_as_state_variable', '', 'act as state variable', & default=self%act_as_state_variable) @@ -365,7 +360,6 @@ subroutine absolute_rate_distributor_initialize(self, configunit) class (type_absolute_rate_distributor), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do/)) call self%register_state_dependency(self%id_target, 'target', '', 'variable to apply sources and sinks to') call self%register_dependency(self%id_w, 'w', '1', 'weights for vertical distribution') call self%register_dependency(self%id_w_int, 'w_int', 'm', 'depth-integrated weights for vertical distribution') @@ -416,7 +410,6 @@ subroutine projected_rate_distributor_initialize(self, configunit) class (type_projected_rate_distributor), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_column/)) call self%register_state_dependency(self%id_target, 'target', '', 'variable to apply sources and sinks to') call self%register_dependency(self%id_h, standard_variables%cell_thickness) call self%register_dependency(self%id_sms, 'sms', '', 'sources-sinks') diff --git a/src/builtin/model_library.F90 b/src/builtin/model_library.F90 index 07f8609d..a1bd4ce9 100644 --- a/src/builtin/model_library.F90 +++ b/src/builtin/model_library.F90 @@ -86,9 +86,9 @@ subroutine create(self,name,model) case ('vertical_integral'); allocate(type_depth_integral::model) case ('bounded_vertical_integral');allocate(type_bounded_depth_integral::model) case ('interior_temporal_mean'); allocate(type_interior_temporal_mean::model) - case ('surface_temporal_mean'); allocate(type_surface_temporal_mean::model) - case ('bottom_temporal_mean'); allocate(type_bottom_temporal_mean::model) - case ('surface_temporal_maximum'); allocate(type_surface_temporal_maximum::model) + case ('surface_temporal_mean'); allocate(type_horizontal_temporal_mean::model) + case ('bottom_temporal_mean'); allocate(type_horizontal_temporal_mean::model) + case ('surface_temporal_maximum'); allocate(type_horizontal_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) ! Add new examples models here @@ -106,7 +106,6 @@ subroutine bottom_field_initialize(self,configunit) class (type_bottom_field), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_bottom/)) call self%register_diagnostic_variable(self%id_result, 'result', '', 'bottom values', source=source_do_bottom) call self%register_dependency(self%id_source, 'source', '', 'interior values') end subroutine bottom_field_initialize @@ -127,7 +126,6 @@ subroutine surface_field_initialize(self,configunit) class (type_surface_field), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_surface/)) call self%register_diagnostic_variable(self%id_result, 'result', '', 'surface values', source=source_do_surface) call self%register_dependency(self%id_source, 'source', '', 'interior values') end subroutine surface_field_initialize @@ -148,7 +146,6 @@ subroutine column_projection_initialize(self,configunit) class (type_column_projection),intent(inout),target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do/)) call self%register_dependency(self%id_source,'source', '', 'horizontal source') call self%register_diagnostic_variable(self%id_result,'result', '', 'interior result') end subroutine column_projection_initialize diff --git a/src/builtin/relaxation.F90 b/src/builtin/relaxation.F90 index 65a8ce43..ec6dd554 100644 --- a/src/builtin/relaxation.F90 +++ b/src/builtin/relaxation.F90 @@ -26,7 +26,6 @@ subroutine interior_relaxation_initialize(self,configunit) class (type_interior_relaxation),intent(inout),target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do/)) call self%register_state_dependency(self%id_original, 'original', '', 'variable that is to be relaxed') call self%register_dependency(self%id_target, 'target', '', 'target to relax towards') call self%get_parameter(self%rate_is_variable, 'rate_is_variable', '', 'use variable relaxation rate', default=.false.) diff --git a/src/builtin/scale.F90 b/src/builtin/scale.F90 index 037d1dc2..95a17b9a 100644 --- a/src/builtin/scale.F90 +++ b/src/builtin/scale.F90 @@ -2,6 +2,7 @@ module fabm_builtin_scale use fabm_types + use fabm_builtin_source, only: copy_fluxes, copy_horizontal_fluxes implicit none @@ -9,24 +10,29 @@ module fabm_builtin_scale public type_scaled_interior_variable, type_scaled_horizontal_variable - type, extends(type_base_model) :: type_scaled_interior_variable + type, extends(type_base_model) :: type_scaled_variable + real(rk) :: weight = 1.0_rk + real(rk) :: offset = 0.0_rk + character(len=attribute_length) :: units = '' + real(rk) :: missing_value = -2.e20_rk + logical :: act_as_state_variable = .false. + logical :: include_background = .false. + integer :: result_output = output_instantaneous + end type + + type, extends(type_scaled_variable) :: type_scaled_interior_variable type (type_dependency_id) :: id_source type (type_diagnostic_variable_id) :: id_result - real(rk) :: weight = 1.0_rk - real(rk) :: offset = 0.0_rk - logical :: include_background = .false. contains procedure :: initialize => scaled_interior_variable_initialize procedure :: do => scaled_interior_variable_do procedure :: after_coupling => scaled_interior_variable_after_coupling end type - type, extends(type_base_model) :: type_scaled_horizontal_variable + type, extends(type_scaled_variable) :: type_scaled_horizontal_variable type (type_horizontal_dependency_id) :: id_source type (type_horizontal_diagnostic_variable_id) :: id_result - real(rk) :: weight = 1.0_rk - real(rk) :: offset = 0.0_rk - logical :: include_background = .false. + integer :: domain = domain_horizontal contains procedure :: initialize => scaled_horizontal_variable_initialize procedure :: do_horizontal => scaled_horizontal_variable_do_horizontal @@ -38,7 +44,11 @@ module fabm_builtin_scale subroutine scaled_interior_variable_initialize(self, configunit) class (type_scaled_interior_variable), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do/)) + + call self%register_dependency(self%id_source, 'source', '', 'source variable') + call self%register_diagnostic_variable(self%id_result, 'result', self%units, 'result', & + missing_value=self%missing_value, output=self%result_output, act_as_state_variable=self%act_as_state_variable) + if (self%act_as_state_variable) call copy_fluxes(self, self%id_result, self%id_source%link%target%name, scale_factor=1.0_rk / self%weight) end subroutine subroutine scaled_interior_variable_do(self, _ARGUMENTS_DO_) @@ -66,7 +76,12 @@ end subroutine scaled_interior_variable_after_coupling subroutine scaled_horizontal_variable_initialize(self, configunit) class (type_scaled_horizontal_variable), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_horizontal/)) + + call self%register_dependency(self%id_source, 'source', '', 'source variable') + call self%register_diagnostic_variable(self%id_result, 'result', self%units, 'result', & + missing_value=self%missing_value, output=self%result_output, act_as_state_variable= & + self%act_as_state_variable, source=source_do_horizontal, domain=self%domain) + if (self%act_as_state_variable) call copy_horizontal_fluxes(self, self%id_result, self%id_source%link%target%name, scale_factor=1.0_rk / self%weight) end subroutine subroutine scaled_horizontal_variable_do_horizontal(self, _ARGUMENTS_HORIZONTAL_) diff --git a/src/builtin/source.F90 b/src/builtin/source.F90 index 79acf4d5..e97764b2 100644 --- a/src/builtin/source.F90 +++ b/src/builtin/source.F90 @@ -69,7 +69,6 @@ subroutine constant_surface_flux_initialize(self, configunit) class (type_constant_surface_flux), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_surface/)) call self%register_state_dependency(self%id_target, 'target', 'UNITS m-3', 'target variable') call self%get_parameter(self%flux, 'flux', 'UNITS m-2 s-1', 'flux (positive for into water)') end subroutine constant_surface_flux_initialize @@ -99,7 +98,6 @@ subroutine external_surface_flux_initialize(self, configunit) class (type_external_surface_flux), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_horizontal/)) call self%add_interior_variable('target', 'UNITS m-3', 'target variable', link=self%target) call self%register_surface_flux(self%target, self%id_target_flux, source=source_do_horizontal) call self%register_dependency(self%id_flux, 'flux', 'UNITS m-2 s-1', 'surface flux') @@ -110,7 +108,6 @@ subroutine external_bottom_flux_initialize(self, configunit) class (type_external_bottom_flux), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_horizontal/)) call self%add_interior_variable('target', 'UNITS m-3', 'target variable', link=self%target) call self%register_bottom_flux(self%target, self%id_target_flux, source=source_do_horizontal) call self%register_dependency(self%id_flux, 'flux', 'UNITS m-2 s-1', 'bottom flux') @@ -121,7 +118,6 @@ subroutine interior_source_initialize(self, configunit) class (type_interior_source), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do/)) call self%add_interior_variable('target', 'UNITS m-3', 'target variable', link=self%target) call self%register_source(self%target, self%id_target_sms) call self%register_dependency(self%id_source, 'source', 'UNITS m-3 s-1', 'source') @@ -144,7 +140,6 @@ subroutine bottom_source_initialize(self, configunit) class (type_bottom_source), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_horizontal/)) call self%add_horizontal_variable('target', 'UNITS m-2', 'target variable', domain=domain_bottom, link=self%target) call self%register_bottom_source(self%target, self%id_target_flux, source=source_do_horizontal) call self%register_dependency(self%id_flux, 'source', 'UNITS m-2 s-1', 'source') @@ -155,7 +150,6 @@ subroutine surface_source_initialize(self, configunit) class (type_surface_source), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_horizontal/)) call self%add_horizontal_variable('target', 'UNITS m-2', 'target variable', domain=domain_bottom, link=self%target) call self%register_surface_source(self%target, self%id_target_flux, source=source_do_horizontal) call self%register_dependency(self%id_flux, 'source', 'UNITS m-2 s-1', 'source') diff --git a/src/builtin/sum.F90 b/src/builtin/sum.F90 index c01fed35..7917e2e5 100644 --- a/src/builtin/sum.F90 +++ b/src/builtin/sum.F90 @@ -4,7 +4,6 @@ module fabm_builtin_sum use fabm_types use fabm_builtin_reduction use fabm_builtin_scale - use fabm_builtin_source implicit none @@ -61,7 +60,7 @@ module fabm_builtin_sum end type type, extends(type_base_sum) :: type_horizontal_weighted_sum - class (type_horizontal_standard_variable), pointer :: aggregate_variable + class (type_horizontal_standard_variable), pointer :: aggregate_variable => null() integer :: domain = domain_horizontal type (type_horizontal_dependency_id), allocatable :: id_terms(:) type (type_horizontal_add_id) :: id_result @@ -169,9 +168,6 @@ subroutine add_component_by_link(self, link, weight, include_background) component => append_component(self, weight, include_background) component%link => link - - ! Temporary: also store the name for use in calls to copy_fluxes from add_to_parent - component%name = link%target%name end subroutine add_component_by_link subroutine after_coupling(self) @@ -179,7 +175,6 @@ subroutine after_coupling(self) type (type_component), pointer :: component real(rk) :: background - integer :: i ! If we are not using the actual sum (because we have 0 or 1 components), skip this routine altogether if (.not. self%active) return @@ -187,10 +182,8 @@ subroutine after_coupling(self) ! At this stage, the background values for all variables (if any) are fixed. We can therefore ! compute background contributions already, and add those to the space- and time-invariant offset. background = 0.0_rk - i = 0 component => self%first do while (associated(component)) - i = i + 1 if (component%include_background) then self%offset = self%offset + component%weight * component%link%target%background_values%value else @@ -317,8 +310,6 @@ subroutine weighted_sum_initialize(self, configunit) class (type_weighted_sum_sms_distributor), pointer :: sms_distributor class (type_scaled_interior_variable), pointer :: scaled_variable - call self%register_implemented_routines((/source_do/)) - n = base_initialize(self) allocate(self%id_terms(n)) @@ -332,36 +323,37 @@ subroutine weighted_sum_initialize(self, configunit) component => component%next end do + ! Allow "do" to be called before merge_components + allocate(self%sources(0)) + if (n == 0) then - ! No components at all - the result is a constant + ! No components - create constant field with offset (typically 0) call self%add_interior_variable('result', self%units, 'result', fill_value=self%offset, missing_value=self%missing_value, & output=self%result_output, source=source_constant, link=self%result_link) return elseif (n == 1) then + ! One component only - set up an alias that will be linked to the actual values 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, & - output=output, act_as_state_variable=self%act_as_state_variable, presence=presence_external_required) + output=output, 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) else - ! One component with scale factor other than 1 (or a user-specified requirement NOT to make a direct link to the source variable) + ! One component with scale factor other than 1 - add a child model to perform the scaling allocate(scaled_variable) - call self%add_child(scaled_variable, '*') - call scaled_variable%register_dependency(scaled_variable%id_source, 'source', '', 'source variable') - call request_coupling_to_component(scaled_variable, scaled_variable%id_source%link, self%first) - call scaled_variable%register_diagnostic_variable(scaled_variable%id_result, 'result', self%units, 'result', & - missing_value=self%missing_value, output=output_none, act_as_state_variable=self%act_as_state_variable) + scaled_variable%units = self%units scaled_variable%weight = self%first%weight - scaled_variable%include_background = self%first%include_background scaled_variable%offset = self%offset + scaled_variable%include_background = self%first%include_background + scaled_variable%act_as_state_variable = self%act_as_state_variable + scaled_variable%missing_value = self%missing_value + scaled_variable%result_output = output_none + call self%add_child(scaled_variable, '*') + call request_coupling_to_component(scaled_variable, scaled_variable%id_source%link, self%first) call self%request_coupling(self%result_link, scaled_variable%id_result%link) - if (self%act_as_state_variable) then - ! This scaled variable acts as a state variable. Create a child model to distribute source terms to the original source variable. - call copy_fluxes(scaled_variable, scaled_variable%id_result, self%first%name, scale_factor=1.0_rk / scaled_variable%weight) - if (associated(self%aggregate_variable)) call scaled_variable%add_to_aggregate_variable(self%aggregate_variable, scaled_variable%id_result) - end if + if (self%act_as_state_variable .and. associated(self%aggregate_variable)) call scaled_variable%add_to_aggregate_variable(self%aggregate_variable, scaled_variable%id_result) end if return end if @@ -369,12 +361,11 @@ subroutine weighted_sum_initialize(self, configunit) self%active = .true. call self%add_interior_variable('result', self%units, 'result', fill_value=0.0_rk, missing_value=self%missing_value, & - output=self%result_output, write_index=self%id_result%sum_index, link=self%id_result%link, source=source_do, & - act_as_state_variable=self%act_as_state_variable) + output=self%result_output, write_index=self%id_result%sum_index, link=self%id_result%link, source=source_do) self%result_link => self%id_result%link - if (self%act_as_state_variable) then - ! NB this does not function yet (hence the act_as_state_variable=.false. above) + if (self%act_as_state_variable .and. .false.) then + ! NB this does not function yet ! Auto-generation of result_sms_tot fails and the do routine of type_weighted_sum_sms_distributor is not yet implemented. ! The sum will act as a state variable. Any source terms will have to be distributed over the individual variables that contribute to the sum. @@ -389,7 +380,7 @@ subroutine weighted_sum_initialize(self, configunit) do i = 1, n write (temp,'(i0)') i call sms_distributor%register_state_dependency(sms_distributor%id_targets(i), 'target' // trim(temp), self%units, 'target ' // trim(temp)) - call sms_distributor%request_coupling(sms_distributor%id_targets(i), trim(component%name)) + call request_coupling_to_component(sms_distributor, sms_distributor%id_targets(i)%link, component) sms_distributor%weights(i) = component%weight component => component%next end do @@ -406,6 +397,7 @@ subroutine weighted_sum_merge_components(self, log_unit) n = base_merge_components(self, log_unit) ! Put all ids and weights in a single array to minimize memory bottlenecks when computing sum + deallocate(self%sources) allocate(self%sources(n)) component => self%first do i = 1, n @@ -441,8 +433,6 @@ subroutine horizontal_weighted_sum_initialize(self, configunit) character(len=10) :: temp class (type_scaled_horizontal_variable), pointer :: scaled_variable - call self%register_implemented_routines((/source_do_horizontal/)) - n = base_initialize(self) allocate(self%id_terms(n)) @@ -456,37 +446,38 @@ subroutine horizontal_weighted_sum_initialize(self, configunit) component => component%next end do + ! Allow "do_horizontal" to be called before merge_components + allocate(self%sources(0)) + if (n == 0) then - ! No components - link to constant field with offset (typically 0) + ! No components - create constant field with offset (typically 0) call self%add_horizontal_variable('result', self%units, 'result', fill_value=self%offset, missing_value=self%missing_value, & output=self%result_output, source=source_constant, link=self%result_link, domain=self%domain) return elseif (n == 1) then - ! One component only. + ! One component only - set up an alias that will be linked to the actual values 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) + output=output, 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) else - ! One component with scale factor other than 1 + ! One component with scale factor other than 1 - add a child model to perform the scaling allocate(scaled_variable) - call self%add_child(scaled_variable, '*') - call scaled_variable%register_dependency(scaled_variable%id_source, 'source', '', 'source variable') - call request_coupling_to_component(scaled_variable, scaled_variable%id_source%link, self%first) - call scaled_variable%register_diagnostic_variable(scaled_variable%id_result, 'result', self%units, 'result', & - missing_value=self%missing_value, output=output_none, act_as_state_variable= & - self%act_as_state_variable, source=source_do_horizontal, domain=self%domain) + scaled_variable%units = self%units scaled_variable%weight = self%first%weight - scaled_variable%include_background = self%first%include_background scaled_variable%offset = self%offset + scaled_variable%domain = self%domain + scaled_variable%include_background = self%first%include_background + scaled_variable%act_as_state_variable = self%act_as_state_variable + scaled_variable%missing_value = self%missing_value + scaled_variable%result_output = output_none + call self%add_child(scaled_variable, '*') + call request_coupling_to_component(scaled_variable, scaled_variable%id_source%link, self%first) call self%request_coupling(self%result_link, scaled_variable%id_result%link) - if (self%act_as_state_variable) then - call copy_horizontal_fluxes(scaled_variable, scaled_variable%id_result, self%first%name, scale_factor=1.0_rk / scaled_variable%weight) - if (associated(self%aggregate_variable)) call scaled_variable%add_to_aggregate_variable(self%aggregate_variable, scaled_variable%id_result) - end if + if (self%act_as_state_variable .and. associated(self%aggregate_variable)) call scaled_variable%add_to_aggregate_variable(self%aggregate_variable, scaled_variable%id_result) end if return end if @@ -508,6 +499,7 @@ subroutine horizontal_weighted_sum_merge_components(self, log_unit) n = base_merge_components(self, log_unit) ! Put all ids and weights in a single array to minimize memory bottlenecks when computing sum + deallocate(self%sources) allocate(self%sources(n)) component => self%first do i = 1, n diff --git a/src/builtin/time_filter.F90 b/src/builtin/time_filter.F90 index 82385221..e7ed1132 100644 --- a/src/builtin/time_filter.F90 +++ b/src/builtin/time_filter.F90 @@ -1,45 +1,74 @@ #include "fabm_driver.h" +#include "fabm_private.h" module fabm_builtin_time_filter use fabm_types - use fabm_expressions, only: temporal_mean, temporal_maximum + use fabm_global_types implicit none private - public type_interior_temporal_mean, type_surface_temporal_mean, type_bottom_temporal_mean, type_surface_temporal_maximum + !public type_interior_temporal_mean, type_surface_temporal_mean, type_bottom_temporal_mean, type_surface_temporal_maximum + public type_interior_temporal_mean, type_horizontal_temporal_mean, type_horizontal_temporal_maximum - type, extends(type_base_model) :: type_interior_temporal_mean - type (type_diagnostic_variable_id) :: id_meanout - type (type_dependency_id) :: id_source, id_meanin - contains - procedure :: initialize => interior_temporal_mean_initialize - procedure :: do => interior_temporal_mean_do - end type + type, extends(type_global_model) :: type_interior_temporal_mean + real(rk) :: window ! Time period to average over (s) + integer :: n + real(rk) :: missing_value = -2.e20_rk + logical :: use_incomplete_result = .false. - type, extends(type_base_model) :: type_surface_temporal_mean - type (type_surface_diagnostic_variable_id) :: id_meanout - type (type_horizontal_dependency_id) :: id_source, id_meanin + type (type_global_interior_dependency_id) :: source + type (type_global_interior_variable_id) :: mean + type (type_global_interior_variable_id), private :: previous_value, last_exact_mean + type (type_global_interior_variable_id), allocatable, private :: history(:) + type (type_global_scalar_variable_id), private :: previous_time, start_time, icurrent contains - procedure :: initialize => surface_temporal_mean_initialize - procedure :: do_surface => surface_temporal_mean_do_surface + procedure :: initialize => interior_temporal_mean_initialize + procedure :: set_data => interior_temporal_mean_set_data + procedure :: update => interior_temporal_mean_update end type - type, extends(type_base_model) :: type_bottom_temporal_mean - type (type_bottom_diagnostic_variable_id) :: id_meanout - type (type_horizontal_dependency_id) :: id_source, id_meanin + type, extends(type_global_model) :: type_horizontal_temporal_mean + real(rk) :: window ! Time period to average over (s) + integer :: n + real(rk) :: missing_value = -2.e20_rk + logical :: use_incomplete_result = .false. + + type (type_global_horizontal_dependency_id) :: source + type (type_global_horizontal_variable_id) :: mean + type (type_global_horizontal_variable_id), private :: previous_value, last_exact_mean + type (type_global_horizontal_variable_id), allocatable, private :: history(:) + type (type_global_scalar_variable_id), private :: previous_time, start_time, icurrent contains - procedure :: initialize => bottom_temporal_mean_initialize - procedure :: do_bottom => bottom_temporal_mean_do_bottom + procedure :: initialize => horizontal_temporal_mean_initialize + procedure :: set_data => horizontal_temporal_mean_set_data + procedure :: update => horizontal_temporal_mean_update end type - - type, extends(type_base_model) :: type_surface_temporal_maximum - type (type_surface_diagnostic_variable_id) :: id_maxout - type (type_horizontal_dependency_id) :: id_source, id_maxin + ! + !type, extends(type_global_model) :: type_bottom_temporal_mean + ! type (type_bottom_diagnostic_variable_id) :: id_output + ! type (type_horizontal_dependency_id) :: id_input + !contains + ! procedure :: initialize => bottom_temporal_mean_initialize + ! procedure :: set_data => bottom_temporal_mean_set_data + ! procedure :: update => bottom_temporal_mean_update + !end type + + type, extends(type_global_model) :: type_horizontal_temporal_maximum + real(rk) :: window ! Time window to compute running maximum over (s) + integer :: n ! Number of bins to use to cover the period + real(rk) :: missing_value = -2.e20_rk ! Missing value to use until the simulation has covered the window size [period] + + type (type_global_horizontal_dependency_id) :: source + type (type_global_horizontal_variable_id) :: maximum + type (type_global_horizontal_variable_id), private :: previous_value + type (type_global_horizontal_variable_id), allocatable, private :: history(:) + type (type_global_scalar_variable_id), private :: previous_time, start_time, current contains - procedure :: initialize => surface_temporal_maximum_initialize - procedure :: do_surface => surface_temporal_maximum_do_surface + procedure :: initialize => horizontal_temporal_maximum_initialize + procedure :: set_data => horizontal_temporal_maximum_set_data + procedure :: update => horizontal_temporal_maximum_update end type contains @@ -48,112 +77,489 @@ subroutine interior_temporal_mean_initialize(self, configunit) class (type_interior_temporal_mean), intent(inout), target :: self integer, intent(in) :: configunit - real(rk) :: window, missing_value - integer :: n - - call self%get_parameter(window, 'window', 's', 'window size') - call self%get_parameter(n, 'n', '', 'number of bins') - call self%get_parameter(missing_value, 'missing_value', '', 'missing value to until the full window size has been covered', default=-2e20_rk) - - call self%register_dependency(self%id_source, 'source', '', 'variable for which to compute running mean') - call self%register_dependency(self%id_meanin, temporal_mean(self%id_source, period=window, resolution=window / n, missing_value=missing_value)) - call self%register_diagnostic_variable(self%id_meanout, 'mean', '', 'running mean') + integer :: ibin + character(len=10) :: strindex + + if (self%user_created) then + call self%get_parameter(self%window, 'window', 's', 'window size') + call self%get_parameter(self%n, 'n', '', 'number of bins') + call self%get_parameter(self%missing_value, 'missing_value', '', 'missing value to until the full window size has been covered', default=-2e20_rk) + end if + + !call self%register_dependency(self%id_input, 'source', '', 'variable for which to compute running mean') + !call self%register_diagnostic_variable(self%id_output, 'mean', '', 'running mean') + call self%add_interior_variable("source", link=self%source%link, presence=presence_external_required) + + allocate(self%history(self%n + 1)) + do ibin = 1, size(self%history) + write (strindex,'(i0)') ibin + call self%add_interior_variable("bin" // trim(strindex), link=self%history(ibin)%link, source=source_global, output=output_none) + self%history(ibin)%link%target%prefill_value = 0.0_rk + self%history(ibin)%link%target%part_of_state = .true. + end do + call self%add_interior_variable("last", link=self%previous_value%link, source=source_global, output=output_none) + call self%add_interior_variable("last_exact_mean", link=self%last_exact_mean%link, source=source_global, output=output_none) + call self%add_interior_variable("mean", link=self%mean%link, source=source_global, output=output_none) + call self%add_scalar_variable("last_time", link=self%previous_time%link, source=source_global, output=output_none) + call self%add_scalar_variable("start_time", link=self%start_time%link, source=source_global, output=output_none) + call self%add_scalar_variable("icurrent", link=self%icurrent%link, source=source_global, output=output_none) + self%last_exact_mean%link%target%prefill_value = 0.0_rk + self%mean%link%target%prefill_value = self%missing_value + self%start_time%link%target%prefill_value = -huge(self%start_time%p) + self%icurrent%link%target%prefill_value = 1.0_rk + self%previous_value%link%target%part_of_state = .true. + self%previous_time%link%target%part_of_state = .true. + self%start_time%link%target%part_of_state = .true. + self%icurrent%link%target%part_of_state = .true. + self%last_exact_mean%link%target%part_of_state = .true. end subroutine - subroutine interior_temporal_mean_do(self, _ARGUMENTS_DO_) - class (type_interior_temporal_mean), intent(in) :: self - _DECLARE_ARGUMENTS_DO_ - - real(rk) :: value - - _LOOP_BEGIN_ - _GET_(self%id_meanin, value) - _SET_DIAGNOSTIC_(self%id_meanout, value) - _LOOP_END_ - end subroutine - - subroutine surface_temporal_mean_initialize(self, configunit) - class (type_surface_temporal_mean), intent(inout), target :: self + subroutine horizontal_temporal_mean_initialize(self, configunit) + class (type_horizontal_temporal_mean), intent(inout), target :: self integer, intent(in) :: configunit - real(rk) :: window, missing_value - integer :: n - - call self%get_parameter(window, 'window', 's', 'window size') - call self%get_parameter(n, 'n', '', 'number of bins') - call self%get_parameter(missing_value, 'missing_value', '', 'missing value to use until the full window size has been covered', default=-2e20_rk) - - call self%register_dependency(self%id_source, 'source', '', 'variable for which to compute running mean') - call self%register_dependency(self%id_meanin, temporal_mean(self%id_source, period=window, resolution=window / n, missing_value=missing_value)) - call self%register_diagnostic_variable(self%id_meanout, 'mean', '', 'running mean') + integer :: ibin + character(len=10) :: strindex + + if (self%user_created) then + call self%get_parameter(self%window, 'window', 's', 'window size') + call self%get_parameter(self%n, 'n', '', 'number of bins') + call self%get_parameter(self%missing_value, 'missing_value', '', 'missing value to use until the full window size has been covered', default=-2e20_rk) + end if + + !call self%register_dependency(self%id_input, 'source', '', 'variable for which to compute running mean') + !call self%register_diagnostic_variable(self%id_output, 'mean', '', 'running mean') + call self%add_horizontal_variable("source", link=self%source%link, presence=presence_external_required) + + allocate(self%history(self%n + 1)) + do ibin = 1, size(self%history) + write (strindex,'(i0)') ibin + call self%add_horizontal_variable("bin" // trim(strindex), link=self%history(ibin)%link, source=source_global, output=output_none) + self%history(ibin)%link%target%prefill_value = 0.0_rk + self%history(ibin)%link%target%part_of_state = .true. + end do + call self%add_horizontal_variable("last", link=self%previous_value%link, source=source_global, output=output_none) + call self%add_horizontal_variable("last_exact_mean", link=self%last_exact_mean%link, source=source_global, output=output_none) + call self%add_horizontal_variable("mean", link=self%mean%link, source=source_global, output=output_none) + call self%add_scalar_variable("last_time", link=self%previous_time%link, source=source_global, output=output_none) + call self%add_scalar_variable("start_time", link=self%start_time%link, source=source_global, output=output_none) + call self%add_scalar_variable("icurrent", link=self%icurrent%link, source=source_global, output=output_none) + self%last_exact_mean%link%target%prefill_value = 0.0_rk + self%mean%link%target%prefill_value = self%missing_value + self%start_time%link%target%prefill_value = -huge(self%start_time%p) + self%icurrent%link%target%prefill_value = 1.0_rk + self%previous_value%link%target%part_of_state = .true. + self%previous_time%link%target%part_of_state = .true. + self%start_time%link%target%part_of_state = .true. + self%icurrent%link%target%part_of_state = .true. + self%last_exact_mean%link%target%part_of_state = .true. end subroutine + ! + !subroutine bottom_temporal_mean_initialize(self, configunit) + ! class (type_bottom_temporal_mean), intent(inout), target :: self + ! integer, intent(in) :: configunit + ! + ! real(rk) :: window, missing_value + ! integer :: n + ! + ! call self%get_parameter(self%window, 'window', 's', 'window size') + ! call self%get_parameter(self%n, 'n', '', 'number of bins') + ! call self%get_parameter(self%missing_value, 'missing_value', '', 'missing value to use until the full window size has been covered', default=-2e20_rk) + ! + ! call self%register_dependency(self%id_input, 'source', '', 'variable for which to compute running mean') + ! call self%register_diagnostic_variable(self%id_output, 'mean', '', 'running mean') + !end subroutine + + subroutine horizontal_temporal_maximum_initialize(self, configunit) + class (type_horizontal_temporal_maximum), intent(inout), target :: self + integer, intent(in) :: configunit - subroutine surface_temporal_mean_do_surface(self, _ARGUMENTS_DO_SURFACE_) - class (type_surface_temporal_mean), intent(in) :: self - _DECLARE_ARGUMENTS_DO_SURFACE_ - - real(rk) :: value - - _SURFACE_LOOP_BEGIN_ - _GET_SURFACE_(self%id_meanin, value) - _SET_SURFACE_DIAGNOSTIC_(self%id_meanout, value) - _SURFACE_LOOP_END_ + integer :: ibin + character(len=10) :: strindex + + if (self%user_created) then + call self%get_parameter(self%window, 'window', 's', 'window size') + call self%get_parameter(self%n, 'n', '', 'number of bins') + call self%get_parameter(self%missing_value, 'missing_value', '', 'missing value to use until the full window size has been covered', default=-2e20_rk) + end if + + !call self%register_dependency(self%id_input, 'source', '', 'variable for which to compute running maximum') + call self%add_horizontal_variable("source", link=self%source%link, presence=presence_external_required) + + allocate(self%history(self%n)) + do ibin = 1, size(self%history) + write (strindex,'(i0)') ibin + call self%add_horizontal_variable("bin" // trim(strindex), link=self%history(ibin)%link, source=source_global, output=output_none) + self%history(ibin)%link%target%prefill_value = -huge(self%history(ibin)%link%target%prefill_value) + self%history(ibin)%link%target%part_of_state = .true. + end do + call self%add_horizontal_variable("last", link=self%previous_value%link, source=source_global, output=output_none) + call self%add_horizontal_variable("max", link=self%maximum%link, source=source_global, output=output_none) + call self%add_scalar_variable("last_time", link=self%previous_time%link, source=source_global, output=output_none) + call self%add_scalar_variable("start_time", link=self%start_time%link, source=source_global, output=output_none) + call self%add_scalar_variable("icurrent", link=self%current%link, source=source_global, output=output_none) + self%maximum%link%target%prefill_value = self%missing_value + self%start_time%link%target%prefill_value = -huge(self%start_time%p) + self%current%link%target%prefill_value = 1.0_rk + self%maximum%link%target%part_of_state = .true. + self%previous_value%link%target%part_of_state = .true. + self%previous_time%link%target%part_of_state = .true. + self%start_time%link%target%part_of_state = .true. + self%current%link%target%part_of_state = self%n > 1 end subroutine - subroutine bottom_temporal_mean_initialize(self, configunit) - class (type_bottom_temporal_mean), intent(inout), target :: self - integer, intent(in) :: configunit - - real(rk) :: window, missing_value - integer :: n - - call self%get_parameter(window, 'window', 's', 'window size') - call self%get_parameter(n, 'n', '', 'number of bins') - call self%get_parameter(missing_value, 'missing_value', '', 'missing value to use until the full window size has been covered', default=-2e20_rk) - - call self%register_dependency(self%id_source, 'source', '', 'variable for which to compute running mean') - call self%register_dependency(self%id_meanin, temporal_mean(self%id_source, period=window, resolution=window / n, missing_value=missing_value)) - call self%register_diagnostic_variable(self%id_meanout, 'mean', '', 'running mean') + subroutine interior_temporal_mean_set_data(self, store, seconds_per_time_unit) + class (type_interior_temporal_mean), intent(inout) :: self + type (type_store), target :: store + real(rke), intent(in) :: seconds_per_time_unit + + integer :: ibin + + self%source%icatalog = self%source%link%target%catalog_index + self%window = self%window / seconds_per_time_unit + do ibin = 1, size(self%history) + self%history(ibin)%p => store%interior(_PREARG_LOCATION_DIMENSIONS_ self%history(ibin)%link%target%store_index) + end do + self%previous_value%p => store%interior(_PREARG_LOCATION_DIMENSIONS_ self%previous_value%link%target%store_index) + self%last_exact_mean%p => store%interior(_PREARG_LOCATION_DIMENSIONS_ self%last_exact_mean%link%target%store_index) + self%mean%p => store%interior(_PREARG_LOCATION_DIMENSIONS_ self%mean%link%target%store_index) + self%previous_time%p => store%scalar(self%previous_time%link%target%store_index) + self%start_time%p => store%scalar(self%start_time%link%target%store_index) + self%icurrent%p => store%scalar(self%icurrent%link%target%store_index) end subroutine - subroutine bottom_temporal_mean_do_bottom(self, _ARGUMENTS_DO_BOTTOM_) - class (type_bottom_temporal_mean), intent(in) :: self - _DECLARE_ARGUMENTS_DO_BOTTOM_ - - real(rk) :: value - - _BOTTOM_LOOP_BEGIN_ - _GET_BOTTOM_(self%id_meanin, value) - _SET_BOTTOM_DIAGNOSTIC_(self%id_meanout, value) - _BOTTOM_LOOP_END_ + subroutine interior_temporal_mean_update(self, catalog _POSTARG_LOCATION_RANGE_, time) + class (type_interior_temporal_mean), intent(in) :: self + type (type_catalog), intent(in) :: catalog + _DECLARE_ARGUMENTS_LOCATION_RANGE_ + real(rke), optional, intent(in) :: time + + real(rke) _ATTRIBUTES_GLOBAL_, pointer :: value + real(rke) :: dt, w, dt_bin, scale, bin_end_time + integer :: icurrent, icurrentbin, ioldest + _DECLARE_LOCATION_ + + value => catalog%interior(self%source%icatalog)%p + if (.not. associated(value)) call self%fatal_error('interior_temporal_mean_update', 'source pointer of not associated.') + + ! Note that all array processing below uses explicit loops in order to respect + ! any limits on the active domain given by the _LOCATION_RANGE_ argument. + + dt_bin = self%window / self%n + + if (self%start_time%p == -huge(self%start_time%p)) then + self%start_time%p = time + self%previous_time%p = time + end if + + icurrent = self%icurrent%p + icurrentbin = mod(icurrent - 1, self%n + 1) + 1 + + do + ioldest = mod(icurrent, self%n + 1) + 1 + bin_end_time = self%start_time%p + dt_bin * icurrent + if (bin_end_time > time) exit + + ! Linearly interpolate to value at end-of-bin time and add that to the current bin + dt = bin_end_time - self%previous_time%p + w = dt / (time - self%previous_time%p) ! weight for current time (leaving 1-w for previous time) + scale = dt / self%window + _BEGIN_GLOBAL_LOOP_ + self%history(icurrentbin)%p _INDEX_LOCATION_ = self%history(icurrentbin)%p _INDEX_LOCATION_ & + + ((1._rke - 0.5_rke * w) * self%previous_value%p _INDEX_LOCATION_ + 0.5_rke * w * value _INDEX_LOCATION_) & + * scale + _END_GLOBAL_LOOP_ + + if (icurrent > self%n) then + ! We already had a complete history (bins covering the full window size). Add the newly full bin, subtract the oldest bin + _BEGIN_GLOBAL_LOOP_ + self%last_exact_mean%p _INDEX_LOCATION_ = self%last_exact_mean%p _INDEX_LOCATION_ & + - self%history(ioldest)%p _INDEX_LOCATION_ + self%history(icurrentbin)%p _INDEX_LOCATION_ + _END_GLOBAL_LOOP_ + else + ! History is incomplete - just add newly filled bin + _BEGIN_GLOBAL_LOOP_ + self%last_exact_mean%p _INDEX_LOCATION_ = self%last_exact_mean%p _INDEX_LOCATION_ & + + self%history(icurrentbin)%p _INDEX_LOCATION_ + _END_GLOBAL_LOOP_ + end if + + ! Update previous time and value to match end of current bin + self%previous_time%p = bin_end_time + _BEGIN_GLOBAL_LOOP_ + self%previous_value%p _INDEX_LOCATION_ = (1._rke - w) * self%previous_value%p _INDEX_LOCATION_ + w * value _INDEX_LOCATION_ + _END_GLOBAL_LOOP_ + + ! Move to next bin: update indices, end time and empty newly current bin + icurrent = icurrent + 1 + icurrentbin = mod(icurrent - 1, self%n + 1) + 1 + self%icurrent%p = icurrent + _BEGIN_GLOBAL_LOOP_ + self%history(icurrentbin)%p _INDEX_LOCATION_ = 0 + _END_GLOBAL_LOOP_ + end do + + ! Compute average of previous and current value, multiply by time difference, pre-divide by window size, and add to current bin. + scale = 0.5_rke * (time - self%previous_time%p) / self%window + _BEGIN_GLOBAL_LOOP_ + self%history(icurrentbin)%p _INDEX_LOCATION_ = self%history(icurrentbin)%p _INDEX_LOCATION_ & + + scale * (self%previous_value%p _INDEX_LOCATION_ + value _INDEX_LOCATION_) + _END_GLOBAL_LOOP_ + + if (icurrent > self%n) then + ! We have a full history covering at least one window size. Update the running mean. + ! The result is an approximation that assumes linear change over the period covered by the oldest bin. + scale = (time - bin_end_time + dt_bin) / dt_bin + _BEGIN_GLOBAL_LOOP_ + self%mean%p _INDEX_LOCATION_ = self%last_exact_mean%p _INDEX_LOCATION_ & + + self%history(icurrentbin)%p _INDEX_LOCATION_ & + - scale * self%history(ioldest)%p _INDEX_LOCATION_ + _END_GLOBAL_LOOP_ + elseif (self%use_incomplete_result) then + if (self%start_time%p == time) then + ! No results just - just the current (first) point in time + _BEGIN_GLOBAL_LOOP_ + self%mean%p _INDEX_LOCATION_ = value _INDEX_LOCATION_ + _END_GLOBAL_LOOP_ + else + ! Use average so far. The integral has been pre-divided by self%window. + ! Undo this and divide instead by time integrated so far + scale = self%window / (time - self%start_time%p) + _BEGIN_GLOBAL_LOOP_ + self%mean%p _INDEX_LOCATION_ = scale * (self%last_exact_mean%p _INDEX_LOCATION_ & + + self%history(icurrentbin)%p _INDEX_LOCATION_) + _END_GLOBAL_LOOP_ + end if + end if + + ! Store current time and value to enable linear interpolation in subsequent call. + self%previous_time%p = time + _BEGIN_GLOBAL_LOOP_ + self%previous_value%p _INDEX_LOCATION_ = value _INDEX_LOCATION_ + _END_GLOBAL_LOOP_ end subroutine - subroutine surface_temporal_maximum_initialize(self, configunit) - class (type_surface_temporal_maximum), intent(inout), target :: self - integer, intent(in) :: configunit - - real(rk) :: window, missing_value - integer :: n - - call self%get_parameter(window, 'window', 's', 'window size') - call self%get_parameter(n, 'n', '', 'number of bins') - call self%get_parameter(missing_value, 'missing_value', '', 'missing value to use until the full window size has been covered', default=-2e20_rk) - - call self%register_dependency(self%id_source, 'source', '', 'variable for which to compute running maximum') - call self%register_dependency(self%id_maxin, temporal_maximum(self%id_source, period=window, resolution=window / n, missing_value=missing_value)) - call self%register_diagnostic_variable(self%id_maxout, 'max', '', 'running maximum') + subroutine horizontal_temporal_mean_set_data(self, store, seconds_per_time_unit) + class (type_horizontal_temporal_mean), intent(inout) :: self + type (type_store), target :: store + real(rke), intent(in) :: seconds_per_time_unit + + integer :: ibin + + self%source%icatalog = self%source%link%target%catalog_index + self%window = self%window / seconds_per_time_unit + do ibin = 1, size(self%history) + self%history(ibin)%p => store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ self%history(ibin)%link%target%store_index) + end do + self%previous_value%p => store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ self%previous_value%link%target%store_index) + self%last_exact_mean%p => store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ self%last_exact_mean%link%target%store_index) + self%mean%p => store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ self%mean%link%target%store_index) + self%previous_time%p => store%scalar(self%previous_time%link%target%store_index) + self%start_time%p => store%scalar(self%start_time%link%target%store_index) + self%icurrent%p => store%scalar(self%icurrent%link%target%store_index) end subroutine - subroutine surface_temporal_maximum_do_surface(self, _ARGUMENTS_DO_SURFACE_) - class (type_surface_temporal_maximum), intent(in) :: self - _DECLARE_ARGUMENTS_DO_SURFACE_ - - real(rk) :: value + subroutine horizontal_temporal_mean_update(self, catalog _POSTARG_LOCATION_RANGE_, time) + class (type_horizontal_temporal_mean), intent(in) :: self + type (type_catalog), intent(in) :: catalog + _DECLARE_ARGUMENTS_LOCATION_RANGE_ + real(rke), optional, intent(in) :: time + + real(rke) _ATTRIBUTES_GLOBAL_HORIZONTAL_, pointer :: value + real(rke) :: dt, w, dt_bin, scale, bin_end_time + integer :: icurrent, icurrentbin, ioldest + _DECLARE_HORIZONTAL_LOCATION_ + + value => catalog%horizontal(self%source%icatalog)%p + if (.not. associated(value)) call self%fatal_error('interior_temporal_mean_update', 'source pointer of not associated.') + + ! Note that all array processing below uses explicit loops in order to respect + ! any limits on the active domain given by the _HORIZONTAL_LOCATION_RANGE_ argument. + + dt_bin = self%window / self%n + + if (self%start_time%p == -huge(self%start_time%p)) then + self%start_time%p = time + self%previous_time%p = time + end if + + icurrent = self%icurrent%p + icurrentbin = mod(icurrent - 1, self%n + 1) + 1 + + do + ioldest = mod(icurrent, self%n + 1) + 1 + bin_end_time = self%start_time%p + dt_bin * icurrent + if (bin_end_time > time) exit + + ! Linearly interpolate to value at end-of-bin time and add that to the current bin + dt = bin_end_time - self%previous_time%p + w = dt / (time - self%previous_time%p) ! weight for current time (leaving 1-w for previous time) + scale = dt / self%window + _BEGIN_OUTER_VERTICAL_LOOP_ + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ = self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ & + + ((1._rke - 0.5_rke * w) * self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ + 0.5_rke * w * value _INDEX_HORIZONTAL_LOCATION_) & + * scale + _END_OUTER_VERTICAL_LOOP_ + + if (icurrent > self%n) then + ! We already had a complete history (bins covering the full window size). Add the newly full bin, subtract the oldest bin + _BEGIN_OUTER_VERTICAL_LOOP_ + self%last_exact_mean%p _INDEX_HORIZONTAL_LOCATION_ = self%last_exact_mean%p _INDEX_HORIZONTAL_LOCATION_ & + - self%history(ioldest)%p _INDEX_HORIZONTAL_LOCATION_ + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + else + ! History is incomplete - just add newly filled bin + _BEGIN_OUTER_VERTICAL_LOOP_ + self%last_exact_mean%p _INDEX_HORIZONTAL_LOCATION_ = self%last_exact_mean%p _INDEX_HORIZONTAL_LOCATION_ & + + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + end if + + ! Update previous time and value to match end of current bin + self%previous_time%p = bin_end_time + _BEGIN_OUTER_VERTICAL_LOOP_ + self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ = (1._rke - w) * self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ + w * value _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + + ! Move to next bin: update indices, end time and empty newly current bin + icurrent = icurrent + 1 + icurrentbin = mod(icurrent - 1, self%n + 1) + 1 + self%icurrent%p = icurrent + _BEGIN_OUTER_VERTICAL_LOOP_ + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ = 0 + _END_OUTER_VERTICAL_LOOP_ + end do + + ! Compute average of previous and current value, multiply by time difference, pre-divide by window size, and add to current bin. + scale = 0.5_rke * (time - self%previous_time%p) / self%window + _BEGIN_OUTER_VERTICAL_LOOP_ + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ = self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ & + + scale * (self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ + value _INDEX_HORIZONTAL_LOCATION_) + _END_OUTER_VERTICAL_LOOP_ + + if (icurrent > self%n) then + ! We have a full history covering at least one window size. Update the running mean. + ! The result is an approximation that assumes linear change over the period covered by the oldest bin. + scale = (time - bin_end_time + dt_bin) / dt_bin + _BEGIN_OUTER_VERTICAL_LOOP_ + self%mean%p _INDEX_HORIZONTAL_LOCATION_ = self%last_exact_mean%p _INDEX_HORIZONTAL_LOCATION_ & + + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ & + - scale * self%history(ioldest)%p _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + elseif (self%use_incomplete_result) then + if (self%start_time%p == time) then + ! No results just - just the current (first) point in time + _BEGIN_OUTER_VERTICAL_LOOP_ + self%mean%p _INDEX_HORIZONTAL_LOCATION_ = value _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + else + ! Use average so far. The integral has been pre-divided by self%window. + ! Undo this and divide instead by time integrated so far + scale = self%window / (time - self%start_time%p) + _BEGIN_OUTER_VERTICAL_LOOP_ + self%mean%p _INDEX_HORIZONTAL_LOCATION_ = scale * (self%last_exact_mean%p _INDEX_HORIZONTAL_LOCATION_ & + + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_) + _END_OUTER_VERTICAL_LOOP_ + end if + end if + + ! Store current time and value to enable linear interpolation in subsequent call. + self%previous_time%p = time + _BEGIN_OUTER_VERTICAL_LOOP_ + self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ = value _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + end subroutine horizontal_temporal_mean_update + + subroutine horizontal_temporal_maximum_set_data(self, store, seconds_per_time_unit) + class (type_horizontal_temporal_maximum), intent(inout) :: self + type (type_store), target :: store + real(rke), intent(in) :: seconds_per_time_unit + + integer :: ibin + + self%source%icatalog = self%source%link%target%catalog_index + self%window = self%window / seconds_per_time_unit + do ibin = 1, size(self%history) + self%history(ibin)%p => store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ self%history(ibin)%link%target%store_index) + end do + self%previous_value%p =>store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ self%previous_value%link%target%store_index) + self%maximum%p => store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ self%maximum%link%target%store_index) + self%previous_time%p => store%scalar(self%previous_time%link%target%store_index) + self%start_time%p => store%scalar(self%start_time%link%target%store_index) + self%current%p => store%scalar(self%current%link%target%store_index) + end subroutine - _SURFACE_LOOP_BEGIN_ - _GET_SURFACE_(self%id_maxin, value) - _SET_SURFACE_DIAGNOSTIC_(self%id_maxout, value) - _SURFACE_LOOP_END_ + subroutine horizontal_temporal_maximum_update(self, catalog _POSTARG_LOCATION_RANGE_, time) + class (type_horizontal_temporal_maximum), intent(in) :: self + type (type_catalog), intent(in) :: catalog + _DECLARE_ARGUMENTS_LOCATION_RANGE_ + real(rke), optional, intent(in) :: time + + real(rke) _ATTRIBUTES_GLOBAL_HORIZONTAL_, pointer :: value + integer :: ibin, icurrent, icurrentbin + real(rke) :: w, bin_end_time + _DECLARE_HORIZONTAL_LOCATION_ + + value => catalog%horizontal(self%source%icatalog)%p + if (.not. associated(value)) call self%fatal_error('interior_temporal_mean_update', 'source pointer of not associated.') + + ! Note that all array processing below uses explicit loops in order to respect + ! any limits on the active domain given by the _HORIZONTAL_LOCATION_RANGE_ argument. + + if (self%start_time%p == -huge(self%start_time%p)) self%start_time%p = time + + icurrent = self%current%p + icurrentbin = mod(icurrent - 1, self%n) + 1 + do + bin_end_time = self%start_time%p + (self%window / self%n) * icurrent + if (bin_end_time > time) exit + + ! Update previous time and value to match end of current bin. For the latter, linearly interpolate to value at end-of-bin time. + w = (bin_end_time - self%previous_time%p) / (time - self%previous_time%p) ! weight for current time (leaving 1-w for previous time) + self%previous_time%p = bin_end_time + _BEGIN_OUTER_VERTICAL_LOOP_ + self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ = (1._rke - w) * self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ + w * value _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + + ! Complete the current bin by taking the maximum over its previous value and the value at end-of-bin time. + _BEGIN_OUTER_VERTICAL_LOOP_ + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ = max(self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_, & + self%previous_value%p _INDEX_HORIZONTAL_LOCATION_) + _END_OUTER_VERTICAL_LOOP_ + + if (icurrent >= self%n) then + ! We have a complete history - compute the maximum over all bins + _BEGIN_OUTER_VERTICAL_LOOP_ + self%maximum%p _INDEX_HORIZONTAL_LOCATION_ = self%history(1)%p _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + do ibin = 2, self%n + _BEGIN_OUTER_VERTICAL_LOOP_ + self%maximum%p _INDEX_HORIZONTAL_LOCATION_ = max(self%maximum%p _INDEX_HORIZONTAL_LOCATION_, & + self%history(ibin)%p _INDEX_HORIZONTAL_LOCATION_) + _END_OUTER_VERTICAL_LOOP_ + end do + end if + + ! Move to next bin: update indices, end time and set maximum of newly current bin to current value (at start of bin) + icurrent = icurrent + 1 + icurrentbin = mod(icurrent -1, self%n) + 1 + self%current%p = icurrent + _BEGIN_OUTER_VERTICAL_LOOP_ + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ = self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ + end do + + ! Update the maximum of the current bin + _BEGIN_OUTER_VERTICAL_LOOP_ + self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_ = max(self%history(icurrentbin)%p _INDEX_HORIZONTAL_LOCATION_, & + value _INDEX_HORIZONTAL_LOCATION_) + _END_OUTER_VERTICAL_LOOP_ + + ! Store current time and value to enable linear interpolation in subsequent call. + self%previous_time%p = time + _BEGIN_OUTER_VERTICAL_LOOP_ + self%previous_value%p _INDEX_HORIZONTAL_LOCATION_ = value _INDEX_HORIZONTAL_LOCATION_ + _END_OUTER_VERTICAL_LOOP_ end subroutine end module diff --git a/src/builtin/tracer.F90 b/src/builtin/tracer.F90 index 6342f710..556fecbe 100644 --- a/src/builtin/tracer.F90 +++ b/src/builtin/tracer.F90 @@ -32,8 +32,6 @@ subroutine initialize(self, configunit) 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) diff --git a/src/c/fabm_c.F90 b/src/c/fabm_c.F90 index c522969a..22039788 100644 --- a/src/c/fabm_c.F90 +++ b/src/c/fabm_c.F90 @@ -30,6 +30,7 @@ module fabm_c integer, parameter :: INTERIOR_DEPENDENCY = 7 integer, parameter :: HORIZONTAL_DEPENDENCY = 8 integer, parameter :: SCALAR_DEPENDENCY = 9 + integer, parameter :: SCALAR_DIAGNOSTIC_VARIABLE = 10 logical, save :: error_occurred = .false. character(len=:), allocatable, save :: error_message @@ -306,11 +307,12 @@ integer(c_int) function model_count(pmodel) bind(c) end function model_count subroutine get_counts(pmodel, nstate_interior, nstate_surface, nstate_bottom, ndiagnostic_interior, ndiagnostic_horizontal, & - ndependencies_interior, ndependencies_horizontal, ndependencies_scalar, nconserved, nparameters, ncouplings) bind(c) + ndiagnostic_scalar, ndependencies_interior, ndependencies_horizontal, ndependencies_scalar, nconserved, nparameters, & + ncouplings) bind(c) !DIR$ ATTRIBUTES DLLEXPORT :: get_counts type (c_ptr), intent(in), value :: pmodel integer(c_int), intent(out) :: nstate_interior, nstate_surface, nstate_bottom - integer(c_int), intent(out) :: ndiagnostic_interior, ndiagnostic_horizontal + integer(c_int), intent(out) :: ndiagnostic_interior, ndiagnostic_horizontal, ndiagnostic_scalar integer(c_int), intent(out) :: ndependencies_interior, ndependencies_horizontal, ndependencies_scalar integer(c_int), intent(out) :: nconserved, nparameters, ncouplings @@ -323,6 +325,7 @@ subroutine get_counts(pmodel, nstate_interior, nstate_surface, nstate_bottom, nd nstate_bottom = size(model%p%bottom_state_variables) ndiagnostic_interior = size(model%p%interior_diagnostic_variables) ndiagnostic_horizontal = size(model%p%horizontal_diagnostic_variables) + ndiagnostic_scalar = size(model%p%scalar_diagnostic_variables) ndependencies_interior = 0 ndependencies_horizontal = 0 @@ -424,6 +427,8 @@ subroutine get_variable_metadata(pmodel, category, index, length, name, units, l variable => model%p%interior_diagnostic_variables(index) case (HORIZONTAL_DIAGNOSTIC_VARIABLE) variable => model%p%horizontal_diagnostic_variables(index) + case (SCALAR_DIAGNOSTIC_VARIABLE) + variable => model%p%scalar_diagnostic_variables(index) case (CONSERVED_QUANTITY) variable => model%p%conserved_quantities(index) end select @@ -450,6 +455,26 @@ subroutine set_variable_save(pmodel, category, index, value) bind(c) end select end subroutine set_variable_save + function get_variable_part_of_state(pmodel, category, index) bind(c) result(part_of_state) + !DIR$ ATTRIBUTES DLLEXPORT :: get_variable_part_of_state + type (c_ptr), intent(in), value :: pmodel + integer(c_int), intent(in), value :: category, index + integer(c_int) :: part_of_state + + type (type_model_wrapper), pointer :: model + + call c_f_pointer(pmodel, model) + + select case (category) + case (INTERIOR_DIAGNOSTIC_VARIABLE) + part_of_state = logical2int(model%p%interior_diagnostic_variables(index)%part_of_state) + case (HORIZONTAL_DIAGNOSTIC_VARIABLE) + part_of_state = logical2int(model%p%horizontal_diagnostic_variables(index)%part_of_state) + case (SCALAR_DIAGNOSTIC_VARIABLE) + part_of_state = logical2int(model%p%scalar_diagnostic_variables(index)%part_of_state) + end select + end function get_variable_part_of_state + function get_variable(pmodel, category, index) bind(c) result(pvariable) !DIR$ ATTRIBUTES DLLEXPORT :: get_variable type (c_ptr), intent(in), value :: pmodel @@ -475,6 +500,8 @@ function get_variable(pmodel, category, index) bind(c) result(pvariable) variable => model%p%interior_diagnostic_variables(index)%original case (HORIZONTAL_DIAGNOSTIC_VARIABLE) variable => model%p%horizontal_diagnostic_variables(index)%original + case (SCALAR_DIAGNOSTIC_VARIABLE) + variable => model%p%scalar_diagnostic_variables(index)%original case (CONSERVED_QUANTITY) variable => model%p%conserved_quantities(index)%target_hz case (INTERIOR_DEPENDENCY, HORIZONTAL_DEPENDENCY, SCALAR_DEPENDENCY) @@ -493,6 +520,8 @@ function get_variable(pmodel, category, index) bind(c) result(pvariable) end if node => node%next end do + case default + stop 'get_variable: Unknown variable category' end select pvariable = c_null_ptr if (associated(variable)) pvariable = c_loc(variable) @@ -1107,6 +1136,22 @@ function get_horizontal_diagnostic_data(pmodel, index) result(ptr) bind(c) if (associated(pvalue)) ptr = c_loc(pvalue) end function get_horizontal_diagnostic_data + function get_scalar_diagnostic_data(pmodel, index) result(ptr) bind(c) + !DIR$ ATTRIBUTES DLLEXPORT :: get_scalar_diagnostic_data + type (c_ptr), intent(in), value :: pmodel + integer(c_int), intent(in), value :: index + type(c_ptr) :: ptr + + real(rke), pointer :: pvalue + + type (type_model_wrapper), pointer :: model + + call c_f_pointer(pmodel, model) + ptr = c_null_ptr + pvalue => model%p%get_scalar_diagnostic_data(index) + if (associated(pvalue)) ptr = c_loc(pvalue) + end function get_scalar_diagnostic_data + function get_standard_variable_data(pmodel, pstandard_variable, horizontal) result(ptr) bind(c) !DIR$ ATTRIBUTES DLLEXPORT :: get_standard_variable_data type (c_ptr), intent(in), value :: pmodel, pstandard_variable diff --git a/src/fabm.F90 b/src/fabm.F90 index 02dd2244..e9fa5492 100644 --- a/src/fabm.F90 +++ b/src/fabm.F90 @@ -21,13 +21,13 @@ module fabm use fabm_expressions use fabm_driver use fabm_properties - use fabm_builtin_depth_integral use fabm_builtin_reduction use fabm_coupling use fabm_job use fabm_schedule use fabm_debug use fabm_work + use fabm_global_types implicit none @@ -164,24 +164,29 @@ module fabm real(rke) :: initial_value = 0.0_rke end type - !> Metadata for an interior diagnostic variable - type, extends(type_fabm_variable) :: type_fabm_interior_diagnostic_variable - class (type_interior_standard_variable), pointer :: standard_variable => null() - + !> Metadata for a diagnostic variable + type, extends(type_fabm_variable) :: type_fabm_diagnostic_variable !> Whether this variable will be included in output and thus needs to be computed. logical :: save = .false. + logical :: part_of_state = .false. + integer :: source end type + !> Metadata for an interior diagnostic variable + type, extends(type_fabm_diagnostic_variable) :: type_fabm_interior_diagnostic_variable + class (type_interior_standard_variable), pointer :: standard_variable => null() + end type + !> Metadata for a horizontal diagnostic variable - type, extends(type_fabm_variable) :: type_fabm_horizontal_diagnostic_variable + type, extends(type_fabm_diagnostic_variable) :: type_fabm_horizontal_diagnostic_variable class (type_horizontal_standard_variable), pointer :: standard_variable => null() + end type - !> Whether this variable will be included in output and thus needs to be computed. - logical :: save = .false. - - integer :: source + !> Metadata for a scalar [0D] diagnostic variable + type, extends(type_fabm_diagnostic_variable) :: type_fabm_scalar_diagnostic_variable + class (type_global_standard_variable), pointer :: standard_variable => null() end type !> Metadata for a conserved quantity @@ -210,6 +215,7 @@ module fabm type (type_fabm_horizontal_state_variable), allocatable, dimension(:) :: bottom_state_variables type (type_fabm_interior_diagnostic_variable), allocatable, dimension(:) :: interior_diagnostic_variables type (type_fabm_horizontal_diagnostic_variable), allocatable, dimension(:) :: horizontal_diagnostic_variables + type (type_fabm_scalar_diagnostic_variable), allocatable, dimension(:) :: scalar_diagnostic_variables type (type_fabm_conserved_quantity), allocatable, dimension(:) :: conserved_quantities !> @} ! --------------------------------------------------------------------------------------------------------------------------- @@ -258,6 +264,8 @@ module fabm type (type_store) :: store type (type_schedules) :: schedules type (type_domain) :: domain + real (rke) :: seconds_per_time_unit = 0.0_rke + real (rke) :: time = 0.0_rke ! --------------------------------------------------------------------------------------------------------------------------- !> @name Memory caches for exchanging information with biogeochemical model instances !> @{ @@ -379,6 +387,7 @@ module fabm procedure :: get_interior_diagnostic_data procedure :: get_horizontal_diagnostic_data + procedure :: get_scalar_diagnostic_data !> @} ! --------------------------------------------------------------------------------------------------------------------------- !> @name Get variable identifiers @@ -546,10 +555,6 @@ subroutine initialize(self) call self%root%add_horizontal_variable('zero_hz', act_as_state_variable=.true., source=source_constant, & missing_value=0.0_rki, output=output_none) - ! Filter out expressions that FABM can handle itself. - ! The remainder, if any, must be handled by the host model. - call filter_expressions(self) - log_unit = -1 if (self%log) then log_unit = get_free_unit() @@ -653,9 +658,6 @@ subroutine set_domain(self _POSTARG_LOCATION_, seconds_per_time_unit) !! between calls to prepare_inputs(). In turn this enables support for built-in time filters such as moving averages. real(rke), optional, intent(in) :: seconds_per_time_unit - class (type_expression), pointer :: expression - real(rke) :: missing_value - if (self%status < status_initialize_done) call fatal_error('set_domain', 'initialize has not yet been called on this model object.') if (self%status >= status_set_domain_done) call fatal_error('set_domain', 'set_domain has already been called on this model object.') self%status = status_set_domain_done @@ -669,57 +671,7 @@ subroutine set_domain(self _POSTARG_LOCATION_, seconds_per_time_unit) self%domain%horizontal_shape(:) = (/_HORIZONTAL_LOCATION_/) #endif - if (present(seconds_per_time_unit)) then - ! Since the host provides information about time, we will support time filters. - ! These includes moving average and moving maximum filters. - expression => self%root%first_expression - do while (associated(expression)) - select type (expression) - class is (type_interior_temporal_mean) - ! Moving average of interior variable - call self%finalize_outputs_job%request_variable(expression%link%target, store=.true.) - expression%in = expression%link%target%catalog_index - expression%period = expression%period / seconds_per_time_unit - allocate(expression%history(_PREARG_LOCATION_ expression%n + 1)) - expression%history = 0.0_rke -#if _FABM_DIMENSION_COUNT_>0 - allocate(expression%previous_value _INDEX_LOCATION_, expression%last_exact_mean _INDEX_LOCATION_, expression%mean _INDEX_LOCATION_) -#endif - expression%last_exact_mean = 0.0_rke - missing_value = expression%missing_value ! To avoid a stack overflow for the next line with ifort 2021.3 - expression%mean = missing_value - call self%link_interior_data(expression%output_name, expression%mean) - class is (type_horizontal_temporal_mean) - ! Moving average of horizontal variable - call self%finalize_outputs_job%request_variable(expression%link%target, store=.true.) - expression%in = expression%link%target%catalog_index - expression%period = expression%period / seconds_per_time_unit - allocate(expression%history(_PREARG_HORIZONTAL_LOCATION_ expression%n + 1)) - expression%history = 0.0_rke -#if _HORIZONTAL_DIMENSION_COUNT_>0 - allocate(expression%previous_value _INDEX_HORIZONTAL_LOCATION_, expression%last_exact_mean _INDEX_HORIZONTAL_LOCATION_, expression%mean _INDEX_HORIZONTAL_LOCATION_) -#endif - expression%last_exact_mean = 0.0_rke - missing_value = expression%missing_value ! To avoid a stack overflow for the next line with ifort 2021.3 - expression%mean = missing_value - call self%link_horizontal_data(expression%output_name, expression%mean) - class is (type_horizontal_temporal_maximum) - ! Moving maximum of horizontal variable - call self%finalize_outputs_job%request_variable(expression%link%target, store=.true.) - expression%in = expression%link%target%catalog_index - expression%period = expression%period / seconds_per_time_unit - allocate(expression%history(_PREARG_HORIZONTAL_LOCATION_ expression%n)) - expression%history = -huge(1.0_rke) -#if _HORIZONTAL_DIMENSION_COUNT_>0 - allocate(expression%previous_value _INDEX_HORIZONTAL_LOCATION_, expression%maximum _INDEX_HORIZONTAL_LOCATION_) -#endif - missing_value = expression%missing_value ! To avoid a stack overflow for the next line with ifort 2021.3 - expression%maximum = missing_value - call self%link_horizontal_data(expression%output_name, expression%maximum) - end select - expression => expression%next - end do - end if + if (present(seconds_per_time_unit)) self%seconds_per_time_unit = seconds_per_time_unit end subroutine set_domain #if _FABM_DIMENSION_COUNT_>0 @@ -941,6 +893,8 @@ subroutine start(self) call cache_create(self%domain, self%cache_fill_values, self%cache_hz) call cache_create(self%domain, self%cache_fill_values, self%cache_vert) + call initialize_global(self%root) + ! For diagnostics that are not needed, set their write index to 0 (rubbish bin) if (self%log) then open(unit=log_unit, file=log_prefix // 'discards.log', action='write', status='replace', iostat=ios) @@ -990,6 +944,16 @@ subroutine start(self) self%status = status_start_done + do ivar = 1, size(self%interior_diagnostic_variables) + self%interior_diagnostic_variables(ivar)%part_of_state = self%interior_diagnostic_variables(ivar)%target%part_of_state + end do + do ivar = 1, size(self%horizontal_diagnostic_variables) + self%horizontal_diagnostic_variables(ivar)%part_of_state = self%horizontal_diagnostic_variables(ivar)%target%part_of_state + end do + do ivar = 1, size(self%scalar_diagnostic_variables) + self%scalar_diagnostic_variables(ivar)%part_of_state = self%scalar_diagnostic_variables(ivar)%target%part_of_state + end do + contains subroutine gather_check_state_data(variables, dat) @@ -1022,6 +986,24 @@ subroutine flag_variables_with_data(variable_list, data_sources) end do end subroutine + recursive subroutine initialize_global(model) + class (type_base_model), intent(inout) :: model + + type (type_model_list_node), pointer :: child + + select type (model) + class is (type_global_model) + call model%set_data(self%store, self%seconds_per_time_unit) + end select + + ! Process children + child => model%children%first + do while (associated(child)) + call initialize_global(child%model) + child => child%next + end do + end subroutine + subroutine report_unfulfilled_dependency(variable) type (type_internal_variable), target :: variable @@ -1606,6 +1588,20 @@ function get_horizontal_diagnostic_data(self, index) result(dat) dat => self%catalog%horizontal(self%horizontal_diagnostic_variables(index)%target%catalog_index)%p end function get_horizontal_diagnostic_data + ! ------------------------------------------------------------------------------------------------------------------------------ + !> Get pointer to data for scalar diagnostic variable + ! ------------------------------------------------------------------------------------------------------------------------------ + function get_scalar_diagnostic_data(self, index) result(dat) + class (type_fabm_model), intent(in) :: self + integer, intent(in) :: index !< variable index + real(rke), pointer :: dat + + _ASSERT_(self%status >= status_start_done, 'get_scalar_diagnostic_data', 'This routine can only be called after model start.') + dat => null() + if (self%scalar_diagnostic_variables(index)%target%catalog_index /= -1) & + dat => self%catalog%scalar(self%scalar_diagnostic_variables(index)%target%catalog_index)%p + end function get_scalar_diagnostic_data + ! ------------------------------------------------------------------------------------------------------------------------------ !> Get pointer to data for interior variable ! ------------------------------------------------------------------------------------------------------------------------------ @@ -2438,6 +2434,8 @@ subroutine process_job(self, job _POSTARG_HORIZONTAL_LOCATION_RANGE_) _VERTICAL_START_ = self%domain%start(_FABM_DEPTH_DIMENSION_INDEX_) _VERTICAL_STOP_ = self%domain%stop(_FABM_DEPTH_DIMENSION_INDEX_) #endif + case (source_global) + call process_global(task, self%catalog _POSTARG_LOCATION_RANGE_, self%time) end select task => task%next end do @@ -2466,7 +2464,6 @@ subroutine prepare_inputs1(self, t) class (type_fabm_model), intent(inout) :: self real(rke), optional, intent(in) :: t - class (type_expression), pointer :: expression _DECLARE_LOCATION_ # if _FABM_DIMENSION_COUNT_ > 0 @@ -2483,27 +2480,8 @@ subroutine prepare_inputs1(self, t) kstop__ = self%domain%stop(3) # endif + if (present(t)) self%time = t call self%process(self%prepare_inputs_job) - - if (present(t)) then - ! The host has provided information about time. Use this to update moving averages, maxima (if any) - expression => self%root%first_expression - do while (associated(expression)) - select type (expression) - class is (type_interior_temporal_mean) - _ASSERT_(associated(self%catalog%interior(expression%in)%p), 'prepare_inputs1', 'source pointer of ' // trim(expression%output_name) // ' not associated.') - call expression%update(t, self%catalog%interior(expression%in)%p _POSTARG_LOCATION_RANGE_) - class is (type_horizontal_temporal_mean) - _ASSERT_(associated(self%catalog%horizontal(expression%in)%p), 'prepare_inputs1', 'source pointer of ' // trim(expression%output_name) // ' not associated.') - call expression%update(t, self%catalog%horizontal(expression%in)%p _POSTARG_HORIZONTAL_LOCATION_RANGE_) - class is (type_horizontal_temporal_maximum) - _ASSERT_(associated(self%catalog%horizontal(expression%in)%p), 'prepare_inputs1', 'source pointer of ' // trim(expression%output_name) // ' not associated.') - call expression%update(t, self%catalog%horizontal(expression%in)%p _POSTARG_HORIZONTAL_LOCATION_RANGE_) - end select - expression => expression%next - end do - end if - end subroutine prepare_inputs1 subroutine prepare_inputs2(self, t, year, month, day, seconds) @@ -2530,9 +2508,10 @@ subroutine classify_variables(self) type (type_fabm_horizontal_state_variable), pointer :: hz_statevar type (type_fabm_interior_diagnostic_variable), pointer :: diagvar type (type_fabm_horizontal_diagnostic_variable), pointer :: hz_diagvar + type (type_fabm_scalar_diagnostic_variable), pointer :: sc_diagvar type (type_fabm_conserved_quantity), pointer :: consvar type (type_internal_variable), pointer :: object - integer :: nstate, nstate_bot, nstate_surf, ndiag, ndiag_hz, ncons + integer :: nstate, nstate_bot, nstate_surf, ndiag, ndiag_hz, ndiag_sc, ncons type (type_aggregate_variable_list) :: aggregate_variable_list type (type_aggregate_variable), pointer :: aggregate_variable @@ -2550,6 +2529,12 @@ subroutine classify_variables(self) link => link%next end do + link => self%root%links%first + do while (associated(link)) + if (link%target%source == source_global) call self%variable_register%add_to_store(link%target) + link => link%next + end do + ! Get list of conserved quantities (map to universal=domain-independent variables where possible) aggregate_variable_list = collect_aggregate_variables(self%root) aggregate_variable => aggregate_variable_list%first @@ -2605,6 +2590,7 @@ subroutine classify_variables(self) nstate_bot = 0 nstate_surf = 0 ndiag_hz = 0 + ndiag_sc = 0 link => self%links_postcoupling%first do while (associated(link)) object => link%target @@ -2626,7 +2612,7 @@ subroutine classify_variables(self) object%source = source_unknown end select elseif (object%source /= source_unknown .or. .not. associated(link%target, link%original)) then - ! Interior diagnostic variable + ! Interior diagnostic variable or coupled variale marked with "output_always_available" ndiag = ndiag + 1 end if case (domain_horizontal, domain_surface, domain_bottom) @@ -2651,9 +2637,14 @@ subroutine classify_variables(self) object%source = source_unknown end select elseif (object%source /= source_unknown .or. .not. associated(link%target, link%original)) then - ! Horizontal diagnostic variable + ! Horizontal diagnostic variable or coupled variale marked with "output_always_available" ndiag_hz = ndiag_hz + 1 end if + case (domain_scalar) + if ((object%source /= source_unknown .and. object%source /= source_state) .or. .not. associated(link%target, link%original)) then + ! Scalar diagnostic variable or coupled variale marked with "output_always_available" + ndiag_sc = ndiag_sc + 1 + end if end select link => link%next end do @@ -2664,6 +2655,7 @@ subroutine classify_variables(self) allocate(self%surface_state_variables (nstate_surf)) allocate(self%interior_diagnostic_variables (ndiag)) allocate(self%horizontal_diagnostic_variables(ndiag_hz)) + allocate(self%scalar_diagnostic_variables (ndiag_sc)) allocate(self%get_interior_sources_job%arg1_sources(nstate)) allocate(self%get_surface_sources_job%arg1_sources(nstate), self%get_surface_sources_job%arg2_sources(nstate_surf)) @@ -2676,6 +2668,7 @@ subroutine classify_variables(self) nstate_bot = 0 nstate_surf = 0 ndiag_hz = 0 + ndiag_sc = 0 link => self%links_postcoupling%first do while (associated(link)) object => link%target @@ -2702,7 +2695,7 @@ subroutine classify_variables(self) call object%movement_sum%target%write_indices%append(self%get_vertical_movement_job%arg1_sources(nstate)) end if elseif (object%source /= source_unknown .or. .not. associated(link%target, link%original)) then - ! Interior diagnostic variable + ! Interior diagnostic variable or coupled variable marked with "output_always_available" ndiag = ndiag + 1 diagvar => self%interior_diagnostic_variables(ndiag) call copy_variable_metadata(link%original, diagvar) @@ -2742,7 +2735,7 @@ subroutine classify_variables(self) hz_statevar%initial_value = object%initial_value end if elseif (object%source /= source_unknown .or. .not. associated(link%target, link%original)) then - ! Horizontal diagnostic variable + ! Horizontal diagnostic variable or coupled variable marked with "output_always_available" ndiag_hz = ndiag_hz + 1 hz_diagvar => self%horizontal_diagnostic_variables(ndiag_hz) call copy_variable_metadata(link%original, hz_diagvar) @@ -2756,6 +2749,22 @@ subroutine classify_variables(self) hz_diagvar%source = object%source hz_diagvar%target => object end if + case (domain_scalar) + if ((object%source /= source_unknown .and. object%source /= source_state) .or. .not. associated(link%target, link%original)) then + ! Scalar diagnostic variable or coupled variable marked with "output_always_available" + ndiag_sc = ndiag_sc + 1 + sc_diagvar => self%scalar_diagnostic_variables(ndiag_sc) + call copy_variable_metadata(link%original, sc_diagvar) + if (associated(object%standard_variables%first)) then + select type (standard_variable => object%standard_variables%first%p) + class is (type_global_standard_variable) + sc_diagvar%standard_variable => standard_variable + end select + end if + sc_diagvar%save = sc_diagvar%output /= output_none + sc_diagvar%source = object%source + sc_diagvar%target => object + end if end select link => link%next end do @@ -2958,6 +2967,7 @@ subroutine create_store(self) ! Collect missing values in array for faster access. These will be used to fill masked parts of outputs. call collect_fill_values(self%variable_register%store%interior, self%store%interior_fill_value, use_missing=.false.) call collect_fill_values(self%variable_register%store%horizontal, self%store%horizontal_fill_value, use_missing=.false.) + call collect_fill_values(self%variable_register%store%scalar, self%store%scalar_fill_value, use_missing=.false.) call collect_fill_values(self%variable_register%store%interior, self%store%interior_missing_value, use_missing=.true.) call collect_fill_values(self%variable_register%store%horizontal, self%store%horizontal_missing_value, use_missing=.true.) @@ -2967,7 +2977,7 @@ subroutine create_store(self) variable_node => self%variable_register%catalog%interior%first do while (associated(variable_node)) if (variable_node%target%store_index > 0) then - ! Note: we first assign to the pointer below to ensure ifort 15 recognizes its contiguity when _FABM_CONTIGUOUS is set + ! Note: we first assign to the pointer below to ensure ifort 15 recognizes its contiguity when _FABM_CONTIGUOUS_ is set pdata => self%store%interior(_PREARG_LOCATION_DIMENSIONS_ variable_node%target%store_index) call self%link_interior_data(variable_node%target, pdata, source=data_source_fabm) end if @@ -2976,12 +2986,19 @@ subroutine create_store(self) variable_node => self%variable_register%catalog%horizontal%first do while (associated(variable_node)) if (variable_node%target%store_index > 0) then - ! Note: we first assign to the pointer below to ensure ifort 15 recognizes its contiguity when _FABM_CONTIGUOUS is set + ! Note: we first assign to the pointer below to ensure ifort 15 recognizes its contiguity when _FABM_CONTIGUOUS_ is set pdata_hz => self%store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ variable_node%target%store_index) call self%link_horizontal_data(variable_node%target, pdata_hz, source=data_source_fabm) end if variable_node => variable_node%next end do + variable_node => self%variable_register%catalog%scalar%first + do while (associated(variable_node)) + if (variable_node%target%store_index > 0) then + call self%link_scalar(variable_node%target, self%store%scalar(variable_node%target%store_index), source=data_source_fabm) + end if + variable_node => variable_node%next + end do contains @@ -2989,6 +3006,7 @@ subroutine allocate_store(_LOCATION_) _DECLARE_ARGUMENTS_LOCATION_ allocate(self%store%interior(_PREARG_LOCATION_ 0:self%variable_register%store%interior%count)) allocate(self%store%horizontal(_PREARG_HORIZONTAL_LOCATION_ 0:self%variable_register%store%horizontal%count)) + allocate(self%store%scalar(0:self%variable_register%store%scalar%count)) end subroutine subroutine collect_fill_values(variable_list, values, use_missing) @@ -3030,6 +3048,9 @@ subroutine reset_store(self) fill_value = self%store%horizontal_fill_value(i) self%store%horizontal(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ i) = fill_value end do + do i = 1, self%variable_register%store%scalar%count + self%store%scalar(i) = self%store%scalar_fill_value(i) + end do end subroutine recursive subroutine merge_indices(model, log_unit) @@ -3051,51 +3072,6 @@ recursive subroutine merge_indices(model, log_unit) end do end subroutine merge_indices - subroutine filter_expressions(self) - class (type_fabm_model),intent(inout) :: self - - class (type_expression), pointer :: current, previous, next - class (type_depth_integral), pointer :: integral - class (type_bounded_depth_integral), pointer :: bounded_integral - logical :: filter - - previous => null() - current => self%root%first_expression - do while (associated(current)) - filter = .false. - select type (current) - class is (type_vertical_integral) - if (current%minimum_depth == 0._rki .and. current%maximum_depth == huge(current%maximum_depth)) then - allocate(integral) - else - allocate(bounded_integral) - bounded_integral%minimum_depth = current%minimum_depth - bounded_integral%maximum_depth = current%maximum_depth - integral => bounded_integral - end if - integral%average = current%average - call self%root%add_child(integral, trim(current%output_name) // '_calculator') - call integral%request_coupling(integral%id_input, current%input_name) - call self%root%request_coupling(current%output_name, integral%id_output%link%target%name) - filter = .true. - end select - - ! If FABM handles this expression internally, remove it from the list. - next => current%next - if (filter) then - if (associated(previous)) then - previous%next => next - else - self%root%first_expression => next - end if - deallocate(current) - else - previous => current - end if - current => next - end do - end subroutine filter_expressions - function get_variable_by_name(self, name, domain) result(variable) class (type_fabm_model), intent(in) :: self character(len=*), intent(in) :: name diff --git a/src/fabm_coupling.F90 b/src/fabm_coupling.F90 index 32259c71..af5903ef 100644 --- a/src/fabm_coupling.F90 +++ b/src/fabm_coupling.F90 @@ -161,7 +161,7 @@ recursive subroutine freeze(self) type (type_model_list_node), pointer :: node - self%frozen = .true. + call self%freeze() node => self%children%first do while (associated(node)) call freeze(node%model) @@ -248,7 +248,7 @@ subroutine collect_user_specified_couplings(self) target_name = self%couplings%get_string(trim(link%name), trim(link%original%long_name), units=trim(link%original%units), default='', display=display) if (len(target_name) >= 4) then if (target_name(len(target_name) - 3:len(target_name)) == '@old') then - call set_dependency_flag(link%original, source=-1, flag=dependency_flag_stale) + call set_dependency_flag(link%original, flag=dependency_flag_stale) target_name = target_name(:len(target_name) - 4) end if end if @@ -661,7 +661,7 @@ recursive subroutine create_aggregate_models(self) 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%act_as_state_variable = associated(aggregate_variable_access%link%target%sms_list%first) sum%result_output = output_none contributing_variable => aggregate_variable%first_contributing_variable diff --git a/src/fabm_expressions.F90 b/src/fabm_expressions.F90 index 39e514ac..bba530db 100644 --- a/src/fabm_expressions.F90 +++ b/src/fabm_expressions.F90 @@ -1,5 +1,4 @@ #include "fabm_driver.h" -#include "fabm_private.h" ! This module define standard expressions that can be used by biogeochemical models. @@ -7,88 +6,55 @@ module fabm_expressions use fabm_types use fabm_driver + use fabm_builtin_depth_integral + use fabm_builtin_time_filter implicit none private public temporal_mean, temporal_maximum, vertical_mean, vertical_integral - public type_interior_temporal_mean, type_horizontal_temporal_mean, type_horizontal_temporal_maximum, type_vertical_integral - type, extends(type_interior_expression) :: type_interior_temporal_mean + type, extends(type_interior_expression) :: type_interior_temporal_mean_expression real(rk) :: period ! Time period to average over (s) integer :: n real(rk) :: missing_value = -2.e20_rk logical :: use_incomplete_result = .false. - type (type_link), pointer :: link => null() - integer :: in = -1 - - real(rk), private :: previous_time, bin_end_time - integer, private :: ioldest = -1 - integer, private :: icurrent = -1 - logical, private :: complete = .false. - real(rke), allocatable _DIMENSION_GLOBAL_PLUS_1_ :: history -#if _FABM_DIMENSION_COUNT_>0 - real(rke), allocatable _DIMENSION_GLOBAL_ :: previous_value, last_exact_mean, mean -#else - real(rke) :: previous_value, last_exact_mean, mean -#endif + type (type_link), pointer :: source => null() contains - procedure :: update => interior_temporal_mean_update + procedure :: resolve => interior_temporal_mean_resolve end type - type, extends(type_horizontal_expression) :: type_horizontal_temporal_mean + type, extends(type_horizontal_expression) :: type_horizontal_temporal_mean_expression real(rk) :: period ! Time period to average over (s) integer :: n real(rk) :: missing_value = 0.0_rk logical :: use_incomplete_result = .false. - type (type_link), pointer :: link => null() - integer :: in = -1 - - real(rk), private :: previous_time, bin_end_time - integer, private :: ioldest = -1 - integer, private :: icurrent = -1 - logical, private :: complete = .false. - real(rke), allocatable _DIMENSION_GLOBAL_HORIZONTAL_PLUS_1_ :: history -#if _HORIZONTAL_DIMENSION_COUNT_>0 - real(rke), allocatable _DIMENSION_GLOBAL_HORIZONTAL_ :: previous_value, last_exact_mean, mean -#else - real(rke) :: previous_value, last_exact_mean, mean -#endif + type (type_link), pointer :: source => null() contains - procedure :: update => horizontal_temporal_mean_update + procedure :: resolve => horizontal_temporal_mean_resolve end type - type, extends(type_horizontal_expression) :: type_horizontal_temporal_maximum + type, extends(type_horizontal_expression) :: type_horizontal_temporal_maximum_expression real(rk) :: period ! Time window to compute running maximum over (s) integer :: n ! Number of bins to use to cover the period real(rk) :: missing_value = -2.e20_rk ! Missing value to use until the simulation has covered the window size [period] - type (type_link), pointer :: link => null() - integer :: in = -1 - - real(rk), private :: previous_time, bin_end_time - integer :: icurrent = -1 - logical, private :: complete = .false. - real(rke), allocatable _DIMENSION_GLOBAL_HORIZONTAL_PLUS_1_ :: history -#if _HORIZONTAL_DIMENSION_COUNT_>0 - real(rke), allocatable _DIMENSION_GLOBAL_HORIZONTAL_ :: previous_value, maximum -#else - real(rke) :: previous_value, maximum -#endif + type (type_link), pointer :: source => null() contains - procedure :: update => horizontal_temporal_maximum_update + procedure :: resolve => horizontal_temporal_maximum_resolve end type type, extends(type_horizontal_expression) :: type_vertical_integral real(rk) :: minimum_depth = 0.0_rk ! Depth below surface in m (positive) real(rk) :: maximum_depth = huge(1.0_rk) ! Depth below surface in m (positive) logical :: average = .false. ! Whether to divide the depth integral by water depth, thus computing the vertical average - character(len=attribute_length) :: input_name = '' - type (type_link), pointer :: link => null() + type (type_link), pointer :: source => null() + contains + procedure :: resolve => vertical_integral_resolve end type interface temporal_mean @@ -100,411 +66,185 @@ module fabm_expressions module procedure horizontal_temporal_maximum end interface - interface vertical_mean - module procedure vertical_dependency_mean - module procedure vertical_state_mean - end interface - - interface vertical_integral - module procedure vertical_dependency_integral - module procedure vertical_state_integral - end interface - contains - function vertical_dependency_mean(input, minimum_depth, maximum_depth) result(expression) - type (type_dependency_id), intent(inout), target :: input - real(rk), optional, intent(in) :: minimum_depth,maximum_depth - type (type_vertical_integral) :: expression - expression = vertical_integral(input,minimum_depth,maximum_depth,average=.true.) - end function - - function vertical_state_mean(input, minimum_depth, maximum_depth) result(expression) - type (type_state_variable_id), intent(inout), target :: input - real(rk), optional, intent(in) :: minimum_depth, maximum_depth - type (type_vertical_integral) :: expression - expression = vertical_integral_generic(input,minimum_depth,maximum_depth,average=.true.) - end function - - function vertical_dependency_integral(input, minimum_depth, maximum_depth, average) result(expression) - type (type_dependency_id), intent(inout),target :: input - real(rk), optional, intent(in) :: minimum_depth, maximum_depth - logical, optional, intent(in) :: average - type (type_vertical_integral) :: expression - expression = vertical_integral_generic(input, minimum_depth, maximum_depth, average) - end function - - function vertical_state_integral(input, minimum_depth, maximum_depth, average) result(expression) - type (type_state_variable_id), intent(inout), target :: input - real(rk), optional, intent(in) :: minimum_depth, maximum_depth - logical, optional, intent(in) :: average - type (type_vertical_integral) :: expression - expression = vertical_integral_generic(input, minimum_depth, maximum_depth, average) - end function - - function vertical_integral_generic(input, minimum_depth, maximum_depth, average) result(expression) - class (type_variable_id), intent(inout), target :: input - real(rk), optional, intent(in) :: minimum_depth, maximum_depth - logical, optional, intent(in) :: average - type (type_vertical_integral) :: expression + function vertical_integral(input, minimum_depth, maximum_depth, average) result(base_expression) + class (type_dependency_id), intent(inout), target :: input + real(rk), optional, intent(in) :: minimum_depth, maximum_depth + logical, optional, intent(in) :: average + class (type_horizontal_expression), pointer :: base_expression - character(len=attribute_length) :: postfix + character(len=attribute_length) :: postfix + type (type_vertical_integral), pointer :: expression - if (.not. associated(input%link)) call fatal_error('fabm_expressions::vertical_mean', & + if (.not. associated(input%link)) call fatal_error('fabm_expressions::vertical_integral', & 'Input variable has not been registered yet.') - expression%input_name = input%link%target%name ! Create a name for the expression - postfix = '' if (present(minimum_depth) .and. present(maximum_depth)) then - if (minimum_depth > maximum_depth) call fatal_error('fabm_expressions::vertical_mean', & + if (minimum_depth > maximum_depth) call fatal_error('fabm_expressions::vertical_integral', & 'Minimum depth exceeds maximum depth.') write (postfix,'(a,i0,a,i0,a)') '_between_', int(minimum_depth), '_m_and_', int(maximum_depth), '_m' elseif (present(minimum_depth)) then write (postfix,'(a,i0,a)') '_below_', int(minimum_depth), '_m' elseif (present(maximum_depth)) then write (postfix,'(a,i0,a)') '_above_', int(maximum_depth), '_m' + else + postfix = '' end if - if (present(average)) expression%average = average + allocate(expression) + if (present(average)) expression%average = average if (expression%average) then expression%output_name = 'vertical_mean_' // trim(input%link%name) // trim(postfix) else expression%output_name = 'integral_of_' // trim(input%link%name) // '_wrt_depth' // trim(postfix) end if - - expression%link => input%link + expression%source => input%link if (present(minimum_depth)) expression%minimum_depth = minimum_depth if (present(maximum_depth)) expression%maximum_depth = maximum_depth + base_expression => expression end function - function interior_temporal_mean(input, period, resolution, missing_value) result(expression) + function vertical_mean(input, minimum_depth, maximum_depth) result(base_expression) + class (type_dependency_id), intent(inout), target :: input + real(rk), optional, intent(in) :: minimum_depth,maximum_depth + class (type_horizontal_expression), pointer :: base_expression + base_expression => vertical_integral(input, minimum_depth, maximum_depth, average=.true.) + end function + + function vertical_integral_resolve(self) result(link) + class (type_vertical_integral), intent(inout) :: self + type (type_link), pointer :: link + + class (type_depth_integral), pointer :: integral + class (type_bounded_depth_integral), pointer :: bounded_integral + + if (self%minimum_depth <= 0._rk .and. self%maximum_depth == huge(self%maximum_depth)) then + allocate(integral) + else + allocate(bounded_integral) + bounded_integral%minimum_depth = self%minimum_depth + bounded_integral%maximum_depth = self%maximum_depth + integral => bounded_integral + end if + integral%average = self%average + call self%link%target%owner%add_child(integral, trim(self%output_name) // '_calculator') + call integral%request_coupling(integral%id_input, self%source) + integral%id_output%link%target%output = output_none + link => integral%id_output%link + end function + + function interior_temporal_mean(input, period, resolution, missing_value) result(base_expression) class (type_dependency_id), intent(inout), target :: input real(rk), intent(in) :: period, resolution real(rk), optional, intent(in) :: missing_value - type (type_interior_temporal_mean) :: expression + class (type_interior_expression), pointer :: base_expression + + type (type_interior_temporal_mean_expression), pointer :: expression if (.not. associated(input%link)) call fatal_error('fabm_expressions::interior_temporal_mean', & 'Input variable has not been registered yet.') + allocate(expression) write (expression%output_name,'(i0,a,a,a,i0,a)') int(period), '_s_mean_', trim(input%link%name), '_at_', int(resolution), '_s_resolution' - expression%link => input%link + expression%source => input%link expression%n = nint(period / resolution) expression%period = period expression%use_incomplete_result = .not. present(missing_value) if (present(missing_value)) expression%missing_value = missing_value + base_expression => expression + end function + + function interior_temporal_mean_resolve(self) result(link) + class (type_interior_temporal_mean_expression), intent(inout) :: self + type (type_link), pointer :: link + + class (type_interior_temporal_mean), pointer :: calculator + + allocate(calculator) + calculator%window = self%period + calculator%n = self%n + calculator%missing_value = self%missing_value + calculator%use_incomplete_result = self%use_incomplete_result + call self%link%target%owner%add_child(calculator, trim(self%output_name) // '_calculator') + call calculator%request_coupling(calculator%source%link, self%source) + calculator%mean%link%target%output = output_none + link => calculator%mean%link end function - function horizontal_temporal_mean(input, period, resolution, missing_value) result(expression) + function horizontal_temporal_mean(input, period, resolution, missing_value) result(base_expression) class (type_horizontal_dependency_id), intent(inout), target :: input real(rk), intent(in) :: period, resolution real(rk), optional, intent(in) :: missing_value - type (type_horizontal_temporal_mean) :: expression + class (type_horizontal_expression), pointer :: base_expression + + type (type_horizontal_temporal_mean_expression), pointer :: expression if (.not. associated(input%link)) call fatal_error('fabm_expressions::horizontal_temporal_mean', & 'Input variable has not been registered yet.') + allocate(expression) write (expression%output_name,'(i0,a,a,a,i0,a)') int(period), '_s_mean_', trim(input%link%name), '_at_', int(resolution), '_s_resolution' - expression%link => input%link + expression%source => input%link expression%n = nint(period / resolution) expression%period = period expression%use_incomplete_result = .not. present(missing_value) if (present(missing_value)) expression%missing_value = missing_value + base_expression => expression + end function + + function horizontal_temporal_mean_resolve(self) result(link) + class (type_horizontal_temporal_mean_expression), intent(inout) :: self + type (type_link), pointer :: link + + class (type_horizontal_temporal_mean), pointer :: calculator + + allocate(calculator) + calculator%window = self%period + calculator%n = self%n + calculator%missing_value = self%missing_value + calculator%use_incomplete_result = self%use_incomplete_result + call self%link%target%owner%add_child(calculator, trim(self%output_name) // '_calculator') + call calculator%request_coupling(calculator%source%link, self%source) + calculator%mean%link%target%output = output_none + link => calculator%mean%link end function - function horizontal_temporal_maximum(input, period, resolution, missing_value) result(expression) + function horizontal_temporal_maximum(input, period, resolution, missing_value) result(base_expression) class (type_horizontal_dependency_id), intent(inout), target :: input real(rk), intent(in) :: period, resolution real(rk), optional, intent(in) :: missing_value - type (type_horizontal_temporal_maximum) :: expression + class (type_horizontal_expression), pointer :: base_expression + + type (type_horizontal_temporal_maximum_expression), pointer :: expression if (.not. associated(input%link)) call fatal_error('fabm_expressions::horizontal_temporal_max', & 'Input variable has not been registered yet.') + allocate(expression) write (expression%output_name,'(i0,a,a,a,i0,a)') int(period), '_s_max_', trim(input%link%name), '_at_', int(resolution), '_s_resolution' - expression%link => input%link + expression%source => input%link expression%n = nint(period / resolution) expression%period = period if (present(missing_value)) expression%missing_value = missing_value + base_expression => expression end function - subroutine interior_temporal_mean_update(self, time, value _POSTARG_LOCATION_RANGE_) - class (type_interior_temporal_mean), intent(inout) :: self - real(rke), intent(in) :: time - real(rke) _ATTRIBUTES_GLOBAL_, intent(in) :: value - _DECLARE_ARGUMENTS_LOCATION_RANGE_ - - real(rke) :: dt, w, dt_bin, scale - _DECLARE_LOCATION_ - - ! Note that all array processing below uses explicit loops in order to respect - ! any limits on the active domain given by the _LOCATION_RANGE_ argument. + function horizontal_temporal_maximum_resolve(self) result(link) + class (type_horizontal_temporal_maximum_expression), intent(inout) :: self + type (type_link), pointer :: link - dt_bin = self%period / self%n + class (type_horizontal_temporal_maximum), pointer :: calculator - if (self%ioldest == -1) then - ! Start of simulation - self%previous_time = time - self%bin_end_time = time + dt_bin - self%icurrent = 1 - self%ioldest = 2 - self%previous_value = 0.0_rke - end if - - do while (time >= self%bin_end_time) - ! Linearly interpolate to value at end-of-bin time and add that to the current bin - dt = self%bin_end_time - self%previous_time - w = dt / (time - self%previous_time) ! weight for current time (leaving 1-w for previous time) - scale = dt / self%period - _BEGIN_GLOBAL_LOOP_ - self%history(_PREARG_LOCATION_ self%icurrent) = self%history(_PREARG_LOCATION_ self%icurrent) & - + ((1._rke - 0.5_rke * w) * self%previous_value _INDEX_LOCATION_ + 0.5_rke * w * value _INDEX_LOCATION_) & - * scale - _END_GLOBAL_LOOP_ - - if (self%complete) then - ! We already had a complete history (bins covering the full window size). Add the newly full bin, subtract the oldest bin - _BEGIN_GLOBAL_LOOP_ - self%last_exact_mean _INDEX_LOCATION_ = self%last_exact_mean _INDEX_LOCATION_ & - - self%history(_PREARG_LOCATION_ self%ioldest) + self%history(_PREARG_LOCATION_ self%icurrent) - _END_GLOBAL_LOOP_ - else - ! History is incomplete - just add newly filled bin - _BEGIN_GLOBAL_LOOP_ - self%last_exact_mean _INDEX_LOCATION_ = self%last_exact_mean _INDEX_LOCATION_ & - + self%history(_PREARG_LOCATION_ self%icurrent) - _END_GLOBAL_LOOP_ - self%complete = self%icurrent == self%n - end if - - ! Update previous time and value to match end of current bin - self%previous_time = self%bin_end_time - _BEGIN_GLOBAL_LOOP_ - self%previous_value _INDEX_LOCATION_ = (1._rke - w) * self%previous_value _INDEX_LOCATION_ + w * value _INDEX_LOCATION_ - _END_GLOBAL_LOOP_ - - ! Move to next bin: update indices, end time and empty newly current bin - self%icurrent = self%ioldest - self%ioldest = self%ioldest + 1 - if (self%ioldest > self%n + 1) self%ioldest = 1 - self%bin_end_time = self%bin_end_time + dt_bin - _BEGIN_GLOBAL_LOOP_ - self%history(_PREARG_LOCATION_ self%icurrent) = 0 - _END_GLOBAL_LOOP_ - end do - - ! Compute average of previous and current value, multiply by time difference, pre-divide by window size, and add to current bin. - scale = 0.5_rke * (time - self%previous_time) / self%period - _BEGIN_GLOBAL_LOOP_ - self%history(_PREARG_LOCATION_ self%icurrent) = self%history(_PREARG_LOCATION_ self%icurrent) & - + scale * (self%previous_value _INDEX_LOCATION_ + value _INDEX_LOCATION_) - _END_GLOBAL_LOOP_ - - if (self%complete) then - ! We have a full history covering at least one window size. Update the running mean. - ! The result is an approximation that assumes linear change over the period covered by the oldest bin. - scale = (time - self%bin_end_time + dt_bin) / dt_bin - _BEGIN_GLOBAL_LOOP_ - self%mean _INDEX_LOCATION_ = self%last_exact_mean _INDEX_LOCATION_ & - + self%history(_PREARG_LOCATION_ self%icurrent) & - - scale * self%history(_PREARG_LOCATION_ self%ioldest) - _END_GLOBAL_LOOP_ - elseif (self%use_incomplete_result) then - if (self%previous_time == time .and. self%icurrent == 1) then - ! No results just - just the current (first) point in time - _BEGIN_GLOBAL_LOOP_ - self%mean _INDEX_LOCATION_ = value _INDEX_LOCATION_ - _END_GLOBAL_LOOP_ - else - ! Use average so far. The integral has been pre-divided by self%period. - ! Undo this and divide instead by time integrated so far - scale = self%period / (self%icurrent * dt_bin + time - self%bin_end_time) - _BEGIN_GLOBAL_LOOP_ - self%mean _INDEX_LOCATION_ = scale * (self%last_exact_mean _INDEX_LOCATION_ & - + self%history(_PREARG_LOCATION_ self%icurrent)) - _END_GLOBAL_LOOP_ - end if - end if - - ! Store current time and value to enable linear interpolation in subsequent call. - self%previous_time = time - _BEGIN_GLOBAL_LOOP_ - self%previous_value _INDEX_LOCATION_ = value _INDEX_LOCATION_ - _END_GLOBAL_LOOP_ - end subroutine - - subroutine horizontal_temporal_mean_update(self, time, value _POSTARG_HORIZONTAL_LOCATION_RANGE_) - class (type_horizontal_temporal_mean), intent(inout) :: self - real(rke), intent(in) :: time - real(rke) _ATTRIBUTES_GLOBAL_HORIZONTAL_, intent(in) :: value - _DECLARE_ARGUMENTS_HORIZONTAL_LOCATION_RANGE_ - - real(rke) :: dt, w, dt_bin, scale - _DECLARE_HORIZONTAL_LOCATION_ - - ! Note that all array processing below uses explicit loops in order to respect - ! any limits on the active domain given by the _HORIZONTAL_LOCATION_RANGE_ argument. - - dt_bin = self%period / self%n - - if (self%ioldest == -1) then - ! Start of simulation - self%previous_time = time - self%bin_end_time = time + dt_bin - self%icurrent = 1 - self%ioldest = 2 - self%previous_value = 0.0_rke - end if - - do while (time >= self%bin_end_time) - ! Linearly interpolate to value at end-of-bin time and add that to the current bin - dt = self%bin_end_time - self%previous_time - w = dt / (time - self%previous_time) ! weight for current time (leaving 1-w for previous time) - scale = dt / self%period - _BEGIN_OUTER_VERTICAL_LOOP_ - self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) = self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) & - + ((1._rke - 0.5_rke * w) * self%previous_value _INDEX_HORIZONTAL_LOCATION_ + 0.5_rke * w * value _INDEX_HORIZONTAL_LOCATION_) & - * scale - _END_OUTER_VERTICAL_LOOP_ - - if (self%complete) then - ! We already had a complete history (bins covering the full window size). Add the newly full bin, subtract the oldest bin - _BEGIN_OUTER_VERTICAL_LOOP_ - self%last_exact_mean _INDEX_HORIZONTAL_LOCATION_ = self%last_exact_mean _INDEX_HORIZONTAL_LOCATION_ & - - self%history(_PREARG_HORIZONTAL_LOCATION_ self%ioldest) + self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) - _END_OUTER_VERTICAL_LOOP_ - else - ! History is incomplete - just add newly filled bin - _BEGIN_OUTER_VERTICAL_LOOP_ - self%last_exact_mean _INDEX_HORIZONTAL_LOCATION_ = self%last_exact_mean _INDEX_HORIZONTAL_LOCATION_ & - + self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) - _END_OUTER_VERTICAL_LOOP_ - self%complete = self%icurrent == self%n - end if - - ! Update previous time and value to match end of current bin - self%previous_time = self%bin_end_time - _BEGIN_OUTER_VERTICAL_LOOP_ - self%previous_value _INDEX_HORIZONTAL_LOCATION_ = (1._rke - w) * self%previous_value _INDEX_HORIZONTAL_LOCATION_ + w * value _INDEX_HORIZONTAL_LOCATION_ - _END_OUTER_VERTICAL_LOOP_ - - ! Move to next bin: update indices, end time and empty newly current bin - self%icurrent = self%ioldest - self%ioldest = self%ioldest + 1 - if (self%ioldest > self%n + 1) self%ioldest = 1 - self%bin_end_time = self%bin_end_time + dt_bin - _BEGIN_OUTER_VERTICAL_LOOP_ - self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) = 0 - _END_OUTER_VERTICAL_LOOP_ - end do - - ! Compute average of previous and current value, multiply by time difference, pre-divide by window size, and add to current bin. - scale = 0.5_rke * (time - self%previous_time) / self%period - _BEGIN_OUTER_VERTICAL_LOOP_ - self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) = self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) & - + scale * (self%previous_value _INDEX_HORIZONTAL_LOCATION_ + value _INDEX_HORIZONTAL_LOCATION_) - _END_OUTER_VERTICAL_LOOP_ - - if (self%complete) then - ! We have a full history covering at least one window size. Update the running mean. - ! The result is an approximation that assumes linear change over the period covered by the oldest bin. - scale = (time - self%bin_end_time + dt_bin) / dt_bin - _BEGIN_OUTER_VERTICAL_LOOP_ - self%mean _INDEX_HORIZONTAL_LOCATION_ = self%last_exact_mean _INDEX_HORIZONTAL_LOCATION_ & - + self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) & - - scale * self%history(_PREARG_HORIZONTAL_LOCATION_ self%ioldest) - _END_OUTER_VERTICAL_LOOP_ - elseif (self%use_incomplete_result) then - if (self%previous_time == time .and. self%icurrent == 1) then - ! No results just - just the current (first) point in time - _BEGIN_OUTER_VERTICAL_LOOP_ - self%mean _INDEX_HORIZONTAL_LOCATION_ = value _INDEX_HORIZONTAL_LOCATION_ - _END_OUTER_VERTICAL_LOOP_ - else - ! Use average so far. The integral has been pre-divided by self%period. - ! Undo this and divide instead by time integrated so far - scale = self%period / (self%icurrent * dt_bin + time - self%bin_end_time) - _BEGIN_OUTER_VERTICAL_LOOP_ - self%mean _INDEX_HORIZONTAL_LOCATION_ = scale * (self%last_exact_mean _INDEX_HORIZONTAL_LOCATION_ & - + self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent)) - _END_OUTER_VERTICAL_LOOP_ - end if - end if - - ! Store current time and value to enable linear interpolation in subsequent call. - self%previous_time = time - _BEGIN_OUTER_VERTICAL_LOOP_ - self%previous_value _INDEX_HORIZONTAL_LOCATION_ = value _INDEX_HORIZONTAL_LOCATION_ - _END_OUTER_VERTICAL_LOOP_ - end subroutine - - subroutine horizontal_temporal_maximum_update(self, time, value _POSTARG_HORIZONTAL_LOCATION_RANGE_) - class (type_horizontal_temporal_maximum), intent(inout) :: self - real(rke), intent(in) :: time - real(rke) _ATTRIBUTES_GLOBAL_HORIZONTAL_, intent(in) :: value - _DECLARE_ARGUMENTS_HORIZONTAL_LOCATION_RANGE_ - - integer :: ibin - real(rke) :: w - _DECLARE_HORIZONTAL_LOCATION_ - - ! Note that all array processing below uses explicit loops in order to respect - ! any limits on the active domain given by the _HORIZONTAL_LOCATION_RANGE_ argument. - - if (self%icurrent == -1) then - ! Start of simulation - self%bin_end_time = time + self%period / self%n - self%icurrent = 1 - end if - - do while (time >= self%bin_end_time) - ! Update previous time and value to match end of current bin. For the latter, linearly interpolate to value at end-of-bin time. - w = (self%bin_end_time - self%previous_time) / (time - self%previous_time) ! weight for current time (leaving 1-w for previous time) - self%previous_time = self%bin_end_time - _BEGIN_OUTER_VERTICAL_LOOP_ - self%previous_value _INDEX_HORIZONTAL_LOCATION_ = (1._rke - w) * self%previous_value _INDEX_HORIZONTAL_LOCATION_ + w * value _INDEX_HORIZONTAL_LOCATION_ - _END_OUTER_VERTICAL_LOOP_ - - ! Complete the current bin by taking the maximum over its previous value and the value at end-of-bin time. - _BEGIN_OUTER_VERTICAL_LOOP_ - self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) = max(self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent), & - self%previous_value _INDEX_HORIZONTAL_LOCATION_) - _END_OUTER_VERTICAL_LOOP_ - - self%complete = self%complete .or. self%icurrent == self%n - if (self%complete) then - ! We have a complete history - compute the maximum over all bins - _BEGIN_OUTER_VERTICAL_LOOP_ - self%maximum _INDEX_HORIZONTAL_LOCATION_ = self%history(_PREARG_HORIZONTAL_LOCATION_ 1) - _END_OUTER_VERTICAL_LOOP_ - do ibin = 2, self%n - _BEGIN_OUTER_VERTICAL_LOOP_ - self%maximum _INDEX_HORIZONTAL_LOCATION_ = max(self%maximum _INDEX_HORIZONTAL_LOCATION_, & - self%history(_PREARG_HORIZONTAL_LOCATION_ ibin)) - _END_OUTER_VERTICAL_LOOP_ - end do - end if - - ! Move to next bin: update indices, end time and set maximum of newly current bin to current value (at start of bin) - self%icurrent = self%icurrent + 1 - if (self%icurrent > self%n) self%icurrent = 1 - self%bin_end_time = self%bin_end_time + self%period / self%n - _BEGIN_OUTER_VERTICAL_LOOP_ - self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) = self%previous_value _INDEX_HORIZONTAL_LOCATION_ - _END_OUTER_VERTICAL_LOOP_ - end do - - ! Update the maximum of the current bin - _BEGIN_OUTER_VERTICAL_LOOP_ - self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent) = max(self%history(_PREARG_HORIZONTAL_LOCATION_ self%icurrent), & - value _INDEX_HORIZONTAL_LOCATION_) - _END_OUTER_VERTICAL_LOOP_ - - ! Store current time and value to enable linear interpolation in subsequent call. - self%previous_time = time - _BEGIN_OUTER_VERTICAL_LOOP_ - self%previous_value _INDEX_HORIZONTAL_LOCATION_ = value _INDEX_HORIZONTAL_LOCATION_ - _END_OUTER_VERTICAL_LOOP_ - end subroutine + allocate(calculator) + calculator%window = self%period + calculator%n = self%n + calculator%missing_value = self%missing_value + call self%link%target%owner%add_child(calculator, trim(self%output_name) // '_calculator') + call calculator%request_coupling(calculator%source%link, self%source) + calculator%maximum%link%target%output = output_none + link => calculator%maximum%link + end function end module fabm_expressions diff --git a/src/fabm_global_types.F90 b/src/fabm_global_types.F90 new file mode 100644 index 00000000..73964529 --- /dev/null +++ b/src/fabm_global_types.F90 @@ -0,0 +1,126 @@ +#include "fabm_driver.h" +#include "fabm_private.h" + +module fabm_global_types + + use fabm_types, only: rke, type_base_model, type_variable_id + + implicit none + + private + + public type_domain, type_catalog, type_store + public type_global_model + public type_global_interior_dependency_id, type_global_horizontal_dependency_id + public type_global_interior_variable_id, type_global_horizontal_variable_id, type_global_scalar_variable_id + + ! -------------------------------------------------------------------------- + ! Derived type for model domain + ! (spatial extent, masks, indices of surface and bottom layers) + ! -------------------------------------------------------------------------- + + type type_domain + ! Information about the model domain + integer :: shape(_FABM_DIMENSION_COUNT_) + integer :: start(_FABM_DIMENSION_COUNT_) + integer :: stop(_FABM_DIMENSION_COUNT_) + integer :: horizontal_shape(_HORIZONTAL_DIMENSION_COUNT_) + +#ifdef _HAS_MASK_ +# ifndef _FABM_HORIZONTAL_MASK_ + _FABM_MASK_TYPE_, pointer _ATTRIBUTES_GLOBAL_ :: mask => null() +# endif + _FABM_MASK_TYPE_, pointer _ATTRIBUTES_GLOBAL_HORIZONTAL_ :: mask_hz => null() +#endif + +#ifdef _FABM_DEPTH_DIMENSION_INDEX_ +# if _FABM_BOTTOM_INDEX_==-1 + integer, pointer _ATTRIBUTES_GLOBAL_HORIZONTAL_ :: bottom_indices => null() +# endif +#endif + end type + + ! -------------------------------------------------------------------------- + ! Derived types for catalog with pointers to all available fields + ! -------------------------------------------------------------------------- + + type type_interior_data_pointer + real(rke), pointer _ATTRIBUTES_GLOBAL_ :: p => null() + end type + + type type_horizontal_data_pointer + real(rke), pointer _ATTRIBUTES_GLOBAL_HORIZONTAL_ :: p => null() + end type + + type type_scalar_data_pointer + real(rke), pointer :: p => null() + end type + + type type_catalog + type (type_interior_data_pointer), allocatable :: interior(:) + type (type_horizontal_data_pointer), allocatable :: horizontal(:) + type (type_scalar_data_pointer), allocatable :: scalar(:) + integer, allocatable :: interior_sources(:) + integer, allocatable :: horizontal_sources(:) + integer, allocatable :: scalar_sources(:) + end type + + ! -------------------------------------------------------------------------- + ! Derived type for variable store + ! (spatially explicit model outputs needed by other BGC modules or host) + ! -------------------------------------------------------------------------- + + type type_store + real(rke), allocatable _DIMENSION_GLOBAL_PLUS_1_ :: interior + real(rke), allocatable _DIMENSION_GLOBAL_HORIZONTAL_PLUS_1_ :: horizontal + real(rke), allocatable, dimension(:) :: scalar + real(rke), allocatable :: interior_fill_value(:) + real(rke), allocatable :: horizontal_fill_value(:) + real(rke), allocatable :: scalar_fill_value(:) + real(rke), allocatable :: interior_missing_value(:) + real(rke), allocatable :: horizontal_missing_value(:) + end type + + + type, extends(type_variable_id) :: type_global_interior_dependency_id + integer :: icatalog = -1 + end type + + type, extends(type_variable_id) :: type_global_horizontal_dependency_id + integer :: icatalog = -1 + end type + + type, extends(type_variable_id) :: type_global_interior_variable_id + real(rke), pointer _ATTRIBUTES_GLOBAL_CONTIGUOUS_ :: p => null() + end type + + type, extends(type_variable_id) :: type_global_horizontal_variable_id + real(rke), pointer _ATTRIBUTES_GLOBAL_HORIZONTAL_CONTIGUOUS_ :: p => null() + end type + + type, extends(type_variable_id) :: type_global_scalar_variable_id + real(rke), pointer :: p => null() + end type + + type, extends(type_base_model) :: type_global_model + contains + procedure :: set_data + procedure :: update + end type + +contains + + subroutine set_data(self, store, seconds_per_time_unit) + class (type_global_model), intent(inout) :: self + type (type_store), target :: store + real(rke), intent(in) :: seconds_per_time_unit + end subroutine + + subroutine update(self, catalog _POSTARG_LOCATION_RANGE_, time) + class (type_global_model), intent(in) :: self + type (type_catalog), intent(in) :: catalog + _DECLARE_ARGUMENTS_LOCATION_RANGE_ + real(rke), optional, intent(in) :: time + end subroutine + +end module diff --git a/src/fabm_graph.F90 b/src/fabm_graph.F90 index ac161c2f..2aea56d3 100644 --- a/src/fabm_graph.F90 +++ b/src/fabm_graph.F90 @@ -394,8 +394,9 @@ recursive function graph_add_call(self, model, source, stack_top) result(node) link => model%links%first do while (associated(link)) - if (index(link%name, '/') == 0 .and. associated(link%original%read_index)) then + if (index(link%name, '/') == 0 .and. (associated(link%original%read_index) .or. (source == source_global .and. link%original%presence == presence_external_required))) then ! This is the model's own variable (not inherited from child model) and the model itself originally requested read access to it. + ! (note: global models/operators do not use read_index, so we detect inputs based on the presence attribute) _ASSERT_(.not. associated(link%target%write_owner), 'graph::add_call', 'BUG: required input variable is co-written.') input_variable => node%inputs%add(link%target) input_variable%update = get_update_flag(link%target, link%original) == dependency_flag_none @@ -471,7 +472,7 @@ recursive subroutine add_call(variable) type (type_output_variable), pointer :: output_variable if (variable%source == source_constant .or. variable%source == source_state .or. variable%source == source_external .or. variable%source == source_unknown) return - _ASSERT_ (.not. variable%write_indices%is_empty(), 'graph_add_variable::add_call', 'Variable "' // trim(variable%name) // '" with source ' // trim(source2string(variable%source)) // ' does not have a write index') + _ASSERT_ (variable%source == source_global .or. .not. variable%write_indices%is_empty(), 'graph_add_variable::add_call', 'Variable "' // trim(variable%name) // '" with source ' // trim(source2string(variable%source)) // ' does not have a write index') node => self%add_call(variable%owner, variable%source, stack_top) output_variable => node%outputs%add(variable) @@ -489,7 +490,7 @@ function source2operation(source) result(operation) integer,intent(in) :: source integer :: operation select case (source) - case (source_do, source_do_column, source_do_bottom, source_do_surface, source_do_horizontal) + case (source_do, source_do_column, source_do_bottom, source_do_surface, source_do_horizontal, source_global) operation = source case (source_get_vertical_movement, source_initialize_state, source_check_state, source_get_light_extinction) operation = source_do diff --git a/src/fabm_job.F90 b/src/fabm_job.F90 index a0c488a9..321b1bdd 100644 --- a/src/fabm_job.F90 +++ b/src/fabm_job.F90 @@ -188,17 +188,6 @@ subroutine task_initialize(self, variable_register, schedules) ! Initialize individual call objects, then collect all input variables in a task-encompassing set. do icall = 1, size(self%calls) - input_variable_node => self%calls(icall)%graph_node%inputs%first - do while (associated(input_variable_node)) - _ASSERT_(.not. input_variable_node%p%target%read_indices%is_empty(), 'task_initialize', 'variable without read indices among inputs') - _ASSERT_(.not. associated(input_variable_node%p%target%write_owner), 'task_initialize', 'write contribution among inputs') - - ! Make sure the variable has an entry in the read register. - call variable_register%add_to_read_cache(input_variable_node%p%target) - - input_variable_node => input_variable_node%next - end do - ! Make sure the pointer to the model has the highest class (and not a base class) ! This is needed because model classes that use inheritance and call base class methods ! may end up with model pointers that are base class specific (and do not reference @@ -218,16 +207,29 @@ subroutine task_initialize(self, variable_register, schedules) end do end if - ! For all output variables that other models are interested in, decide whether to copy their value - ! from the write to read cache [if the other model will be called as part of the same task], - ! of to save it to the persistent data store. - variable_node => self%calls(icall)%graph_node%outputs%first - do while (associated(variable_node)) - call variable_register%add_to_write_cache(variable_node%p%target) - if (variable_node%p%copy_to_cache) call variable_register%add_to_read_cache(variable_node%p%target) - if (variable_node%p%copy_to_store) call variable_register%add_to_store(variable_node%p%target) - variable_node => variable_node%next - end do + if (self%calls(icall)%source /= source_global) then + input_variable_node => self%calls(icall)%graph_node%inputs%first + do while (associated(input_variable_node)) + _ASSERT_(.not. input_variable_node%p%target%read_indices%is_empty(), 'task_initialize', 'variable without read indices among inputs') + _ASSERT_(.not. associated(input_variable_node%p%target%write_owner), 'task_initialize', 'write contribution among inputs') + + ! Make sure the variable has an entry in the read register. + call variable_register%add_to_read_cache(input_variable_node%p%target) + + input_variable_node => input_variable_node%next + end do + + ! For all output variables that other models are interested in, decide whether to copy their value + ! from the write to read cache [if the other model will be called as part of the same task], + ! of to save it to the persistent data store. + variable_node => self%calls(icall)%graph_node%outputs%first + do while (associated(variable_node)) + call variable_register%add_to_write_cache(variable_node%p%target) + if (variable_node%p%copy_to_cache) call variable_register%add_to_read_cache(variable_node%p%target) + if (variable_node%p%copy_to_store) call variable_register%add_to_store(variable_node%p%target) + variable_node => variable_node%next + end do + end if call schedules%attach(self%calls(icall)%model, self%calls(icall)%source, self%calls(icall)%active) end do @@ -412,7 +414,7 @@ subroutine create_persistent_store_commands(commands, domain) do icall = 1, size(self%calls) variable_node => self%calls(icall)%graph_node%outputs%first do while (associated(variable_node)) - if (variable_node%p%copy_to_store .and. iand(variable_node%p%target%domain, domain) /= 0 .and. self%calls(icall)%graph_node%source /= source_constant) then + if (variable_node%p%copy_to_store .and. iand(variable_node%p%target%domain, domain) /= 0 .and. self%calls(icall)%graph_node%source /= source_constant .and. self%calls(icall)%graph_node%source /= source_global) then _ASSERT_(variable_node%p%target%write_indices%value > 0, 'create_prefill_commands', 'Variable ' // trim(variable_node%p%target%name) // ' has copy_to_store set, but it does not have a write cache index.') _ASSERT_(variable_node%p%target%store_index /= store_index_none, 'create_persistent_store_commands', 'Variable ' // trim(variable_node%p%target%name) // ' has copy_to_store set, but it does not have a persistent storage index.') ilast = max(ilast,variable_node%p%target%store_index) @@ -429,7 +431,7 @@ subroutine create_persistent_store_commands(commands, domain) do icall = 1, size(self%calls) variable_node => self%calls(icall)%graph_node%outputs%first do while (associated(variable_node)) - if (variable_node%p%copy_to_store .and. iand(variable_node%p%target%domain,domain) /= 0 .and. self%calls(icall)%graph_node%source /= source_constant) & + if (variable_node%p%copy_to_store .and. iand(variable_node%p%target%domain,domain) /= 0 .and. self%calls(icall)%graph_node%source /= source_constant .and. self%calls(icall)%graph_node%source /= source_global) & commands(variable_node%p%target%store_index) = variable_node%p%target%write_indices%value variable_node => variable_node%next end do @@ -1138,12 +1140,12 @@ subroutine prepare_task(task) do while (associated(input_variable)) final_output = link_cowritten_outputs(input_variable%p%sources) if (associated(final_output%task, task) .and. final_output%icall < icall) then - ! The call that is responsible for computing this input is part of the same task. + ! The call that is responsible for computing this input happens earlier as part of the same task. ! Therefore the output needs to be copied to the read cache. - final_output%variable%copy_to_cache = .true. + if (task%operation /= source_global) final_output%variable%copy_to_cache = .true. else if (associated(final_output%variable)) final_output%variable%copy_to_store = .true. - if (input_variable%p%target%source /= source_constant .and. input_variable%p%target%source /= source_unknown) & + if (input_variable%p%target%source /= source_constant .and. input_variable%p%target%source /= source_unknown .and. task%operation /= source_global) & call task%read_cache_preload%add(input_variable%p%target) end if @@ -1445,7 +1447,11 @@ subroutine job_manager_initialize(self, variable_register, schedules, log_unit, input_variable => graph_node%p%inputs%first do while (associated(input_variable)) if (.not. input_variable%p%update) then + ! Force calculation of this variable (if already calculated here, or earlier, this does nothing) call finalize_job%graph%add_variable(input_variable%p%target, input_variable%p%sources) + + ! Ensure the variable is added to restarts + input_variable%p%target%part_of_state = .true. end if input_variable => input_variable%next end do diff --git a/src/fabm_library.F90.in b/src/fabm_library.F90.in index 754a237a..ec2c9c3e 100644 --- a/src/fabm_library.F90.in +++ b/src/fabm_library.F90.in @@ -40,7 +40,7 @@ contains model => null() - call self%type_base_model_factory%create(name,model) + call self%type_base_model_factory%create(name, model) ! Store name that was used to create this model, so we can re-create it in the future. if (associated(model)) model%type_name = name diff --git a/src/fabm_particle.F90 b/src/fabm_particle.F90 index 12f6b292..132436e5 100644 --- a/src/fabm_particle.F90 +++ b/src/fabm_particle.F90 @@ -55,7 +55,6 @@ module fabm_particle type, extends(type_coupling_task) :: type_coupling_from_model type (type_model_reference), pointer :: model_reference => null() - integer :: access = access_read class (type_particle_model), pointer :: owner => null() contains procedure :: resolve => coupling_from_model_resolve @@ -240,8 +239,6 @@ subroutine request_coupling_to_model_generic(self, link, target_model, target_mo class is (type_domain_specific_standard_variable) coupling%target_standard_variable => standard_variable end select - if (coupling%link%target%source == source_state .or. coupling%link%target%fake_state_variable) & - coupling%access = access_state end if if (present(target_model)) then coupling%model_reference => target_model%reference @@ -402,7 +399,7 @@ recursive subroutine before_coupling(self) type (type_model_reference), pointer :: reference class (type_base_model), pointer :: model - ! For model references that include the model state, buid the corresponding lists + ! For model references that include the model state, build the corresponding lists ! of state variable identifiers now, and request their coupling to the referenced model. ! This must be done before the actual [variable] coupling starts, and therefore has ! to happen here, in before_coupling. @@ -434,7 +431,7 @@ function coupling_from_model_resolve(self) result(link) link => null() elseif (associated(self%target_standard_variable)) then ! Coupling to a standard [aggregate] variable - link => get_aggregate_variable_access(model, self%target_standard_variable, self%access) + link => get_aggregate_variable_access(model, self%target_standard_variable) else ! Coupling to a named variable link => model%find_link(trim(self%target_name)) diff --git a/src/fabm_types.F90 b/src/fabm_types.F90 index 1197ff73..78d2e701 100644 --- a/src/fabm_types.F90 +++ b/src/fabm_types.F90 @@ -103,7 +103,8 @@ module fabm_types source_get_drag = 15, & source_get_albedo = 16, & source_external = 17, & - source_state = 18 + source_state = 18, & + source_global = 19 integer, parameter, public :: presence_internal = 1, & presence_external_required = 2, & @@ -113,11 +114,6 @@ module fabm_types prefill_constant = -1, & prefill_previous_value = -2 - integer, parameter, public :: access_none = 0, & - access_read = 1, & - access_add_source = 2, & - access_state = ior(access_read, access_add_source) - integer, parameter, public :: store_index_none = -1 integer, parameter, public :: operator_assign = 0, & @@ -376,6 +372,7 @@ module fabm_types type (type_dependency_flag), allocatable :: dependency_flags(:) logical :: fake_state_variable = .false. + logical :: part_of_state = .false. ! Only used for interior state variables: logical :: no_precipitation_dilution = .false. @@ -415,18 +412,14 @@ module fabm_types ! variables). ! -------------------------------------------------------------------------- - type, abstract :: type_expression - class (type_expression), pointer :: next => null() - character(len=attribute_length) :: output_name = '' - integer, pointer :: out => null() + type, extends(type_coupling_task) :: type_expression + character(len=attribute_length) :: output_name = '' end type type, abstract, extends(type_expression) :: type_interior_expression - !type (type_interior_data_pointer), pointer :: out => null() end type type, abstract, extends(type_expression) :: type_horizontal_expression - !type (type_horizontal_data_pointer), pointer :: out => null() end type ! -------------------------------------------------------------------------- @@ -480,8 +473,6 @@ module fabm_types type (type_fabm_settings) :: parameters type (type_fabm_settings) :: initialization - class (type_expression), pointer :: first_expression => null() - type (type_coupling_task_list) :: coupling_task_list real(rk) :: dt = 1.0_rk @@ -492,8 +483,6 @@ module fabm_types type (type_add_id) :: extinction_id type (type_horizontal_add_id) :: albedo_id type (type_horizontal_add_id) :: surface_drag_id - - integer, allocatable :: implemented(:) contains ! Procedure for adding child models [during initialization only] @@ -538,7 +527,8 @@ module fabm_types procedure :: set_variable_property_real procedure :: set_variable_property_integer procedure :: set_variable_property_logical - generic :: set_variable_property => set_variable_property_real,set_variable_property_integer,set_variable_property_logical + generic :: set_variable_property => set_variable_property_real, set_variable_property_integer, & + set_variable_property_logical procedure :: add_variable_to_aggregate_variable procedure :: add_constant_to_aggregate_variable @@ -660,7 +650,7 @@ module fabm_types procedure :: after_coupling => base_after_coupling procedure :: implements - procedure :: register_implemented_routines + procedure :: freeze procedure :: finalize => base_finalize @@ -669,6 +659,9 @@ module fabm_types procedure :: get_light_extinction => base_get_light_extinction procedure :: get_drag => base_get_drag procedure :: get_albedo => base_get_albedo + + ! Deprecated as of FABM 3.0 + procedure :: register_implemented_routines end type type_base_model ! ==================================================================================================== @@ -677,7 +670,7 @@ module fabm_types type type_cache ! Number of active items in a single cache line [first dimension of any spatially explicit caches below] - integer :: n = 1 + integer :: n = 0 ! Read cache (separate interior, horizontal, scalar fields). real(rk), allocatable _DIMENSION_SLICE_PLUS_1_ :: read @@ -694,6 +687,7 @@ module fabm_types logical :: valid logical :: set_interior logical :: set_horizontal + logical :: implemented = .true. end type type, extends(type_cache) :: type_interior_cache @@ -759,17 +753,20 @@ subroutine base_initialize(self, configunit) subroutine base_initialize_state(self, _ARGUMENTS_INITIALIZE_STATE_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_INITIALIZE_STATE_ + cache%implemented = .false. end subroutine subroutine base_initialize_horizontal_state(self, _ARGUMENTS_INITIALIZE_HORIZONTAL_STATE_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_INITIALIZE_HORIZONTAL_STATE_ + cache%implemented = .false. end subroutine ! Providing process rates and diagnostics subroutine base_do(self, _ARGUMENTS_DO_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_DO_ + cache%implemented = .false. end subroutine subroutine base_do_ppdd(self, _ARGUMENTS_DO_PPDD_) @@ -781,47 +778,56 @@ subroutine base_do_ppdd(self, _ARGUMENTS_DO_PPDD_) subroutine base_do_bottom(self, _ARGUMENTS_DO_BOTTOM_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_DO_BOTTOM_ + cache%implemented = .false. end subroutine subroutine base_do_bottom_ppdd(self, _ARGUMENTS_DO_BOTTOM_PPDD_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_DO_BOTTOM_PPDD_ + cache%implemented = .false. end subroutine subroutine base_do_surface(self, _ARGUMENTS_DO_SURFACE_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_DO_SURFACE_ + cache%implemented = .false. end subroutine subroutine base_do_horizontal(self, _ARGUMENTS_HORIZONTAL_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_HORIZONTAL_ + cache%implemented = .false. end subroutine subroutine base_do_column(self, _ARGUMENTS_DO_COLUMN_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_DO_COLUMN_ call self%get_light(_ARGUMENTS_DO_COLUMN_) + cache%implemented = .false. end subroutine subroutine base_get_vertical_movement(self, _ARGUMENTS_GET_VERTICAL_MOVEMENT_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_GET_VERTICAL_MOVEMENT_ + cache%implemented = .false. end subroutine subroutine base_check_state(self, _ARGUMENTS_CHECK_STATE_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_CHECK_STATE_ + cache%implemented = .false. end subroutine subroutine base_check_surface_state(self, _ARGUMENTS_CHECK_SURFACE_STATE_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_CHECK_SURFACE_STATE_ + cache%implemented = .false. end subroutine subroutine base_check_bottom_state(self, _ARGUMENTS_CHECK_BOTTOM_STATE_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_CHECK_BOTTOM_STATE_ + cache%implemented = .false. end subroutine recursive subroutine base_finalize(self) @@ -830,7 +836,6 @@ recursive subroutine base_finalize(self) type (type_model_list_node), pointer :: node class (type_base_model), pointer :: child type (type_aggregate_variable_access), pointer :: aggregate_variable_access, next_aggregate_variable_access - class (type_expression), pointer :: expression, next_expression type (type_link), pointer :: link node => self%children%first @@ -851,14 +856,6 @@ recursive subroutine base_finalize(self) end do self%first_aggregate_variable_access => null() - expression => self%first_expression - do while (associated(expression)) - next_expression => expression%next - deallocate(expression) - expression => next_expression - end do - self%first_expression => null() - link => self%links%first do while (associated(link)) if (index(link%name, '/') == 0) then @@ -904,21 +901,25 @@ subroutine finalize_variable(variable) subroutine base_get_light_extinction(self, _ARGUMENTS_GET_EXTINCTION_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_GET_EXTINCTION_ + cache%implemented = .false. end subroutine subroutine base_get_drag(self, _ARGUMENTS_GET_DRAG_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_GET_DRAG_ + cache%implemented = .false. end subroutine subroutine base_get_albedo(self, _ARGUMENTS_GET_ALBEDO_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_GET_ALBEDO_ + cache%implemented = .false. end subroutine subroutine base_get_light(self, _ARGUMENTS_DO_COLUMN_) class (type_base_model), intent(in) :: self _DECLARE_ARGUMENTS_DO_COLUMN_ + cache%implemented = .false. end subroutine function base_get_path(self) result(path) @@ -968,27 +969,62 @@ function implements(self, source) result(is_implemented) integer, intent(in) :: source logical :: is_implemented - integer :: i + integer :: i + type (type_interior_cache) :: interior_cache + type (type_horizontal_cache) :: horizontal_cache - is_implemented = .true. - if (allocated(self%implemented)) then - do i = 1, size(self%implemented) - if (self%implemented(i) == source) return - end do - is_implemented = .false. - end if + allocate(interior_cache%read_scalar(-1:-1), horizontal_cache%read_scalar(-1:-1)) + interior_cache%read_scalar(:) = 0.0_rk + horizontal_cache%read_scalar(:) = 0.0_rk + select case (source) + case (source_initialize_state) + call self%initialize_state(interior_cache) + is_implemented = interior_cache%implemented + case (source_initialize_bottom_state) + call self%initialize_bottom_state(horizontal_cache) + is_implemented = horizontal_cache%implemented + case (source_initialize_surface_state) + call self%initialize_surface_state(horizontal_cache) + is_implemented = horizontal_cache%implemented + case (source_check_state) + call self%check_state(interior_cache) + is_implemented = interior_cache%implemented + case (source_check_bottom_state) + call self%check_bottom_state(horizontal_cache) + is_implemented = horizontal_cache%implemented + case (source_check_surface_state) + call self%check_surface_state(horizontal_cache) + is_implemented = horizontal_cache%implemented + case (source_do) + call self%do(interior_cache) + is_implemented = interior_cache%implemented + case (source_do_surface) + call self%do_surface(horizontal_cache) + is_implemented = horizontal_cache%implemented + case (source_do_bottom) + call self%do_bottom(horizontal_cache) + is_implemented = horizontal_cache%implemented + case (source_get_vertical_movement) + call self%get_vertical_movement(interior_cache) + is_implemented = interior_cache%implemented + case (source_get_light_extinction) + call self%get_light_extinction(interior_cache) + is_implemented = interior_cache%implemented + case (source_get_drag) + call self%get_drag(horizontal_cache) + is_implemented = horizontal_cache%implemented + case (source_get_albedo) + call self%get_albedo(horizontal_cache) + is_implemented = horizontal_cache%implemented + case default + is_implemented = .true. + end select end function subroutine register_implemented_routines(self, sources) class (type_base_model), intent(inout) :: self integer, optional, intent(in) :: sources(:) - if (allocated(self%implemented)) deallocate(self%implemented) - if (present(sources)) then - allocate(self%implemented(size(sources))) - self%implemented(:) = sources - else - allocate(self%implemented(0)) - end if + !call self%log_message('Warning: register_implemented_routines is superfluous as of FABM 3.0 and will be removed in future versions.') end subroutine recursive subroutine add_child(self, model, name, long_name, configunit) @@ -1090,6 +1126,33 @@ recursive subroutine add_child(self, model, name, long_name, configunit) end if end subroutine add_child + subroutine freeze(self) + class (type_base_model), intent(inout) :: self + + type (type_link), pointer :: link + + link => self%links%first + do while (associated(link)) + if (index(link%name, '/') == 0) then + if (link%original%source /= source_unknown .and. link%original%source /= source_state .and. link%original%source /= source_constant .and. link%original%source /= source_do_column) then + if (.not. self%implements(link%original%source)) then + if (link%original%write_operator == operator_add) then + ! Quietly change to no-op - the base class would just not have any effect + link%original%source = source_constant + else + ! Throw an error - the variable will not be given a value + call self%fatal_error('freeze', trim(link%name) // ' is registered with source ' & + // trim(source2string(link%original%source)) // ', but the routine for this is not implemented.') + end if + end if + end if + end if + link => link%next + end do + + self%frozen = .true. + end subroutine + subroutine set_variable_property_real(self, variable, name, value) class (type_base_model), intent(inout) :: self class (type_variable_id), intent(inout) :: variable @@ -1422,12 +1485,12 @@ subroutine request_coupling_lt(self, link, task) if (associated(current_link, link)) exit current_link => current_link%next end do - if (.not.associated(current_link)) call self%fatal_error('request_coupling_ln', & + if (.not. associated(current_link)) call self%fatal_error('request_coupling', & 'Couplings can only be requested for variables that you own yourself.') ! Make sure that the link also points to a variable that we registered ourselves, ! rather than one registered by a child model. - if (index(link%name, '/') /= 0) call self%fatal_error('request_coupling_ln', & + if (index(link%name, '/') /= 0) call self%fatal_error('request_coupling', & 'Couplings can only be requested for variables that you registered yourself, & ¬ inherited ones such as the current ' // trim(link%name) // '.') @@ -1558,7 +1621,7 @@ subroutine integer_pointer_set_extend(self, other) integer :: i if (allocated(other%pointers)) then - do i=1,size(other%pointers) + do i = 1, size(other%pointers) call self%append(other%pointers(i)%p) end do end if @@ -1578,7 +1641,7 @@ subroutine integer_pointer_set_set_value(self, value) integer :: i if (allocated(self%pointers)) then - do i=1,size(self%pointers) + do i = 1, size(self%pointers) self%pointers(i)%p = value end do end if @@ -1619,7 +1682,7 @@ subroutine real_pointer_set_extend(self, other) integer :: i if (allocated(other%pointers)) then - do i=1,size(other%pointers) + do i = 1, size(other%pointers) call self%append(other%pointers(i)%p) end do end if @@ -1632,7 +1695,7 @@ subroutine real_pointer_set_set_value(self, value) integer :: i if (allocated(self%pointers)) then - do i=1,size(self%pointers) + do i = 1, size(self%pointers) self%pointers(i)%p = value end do end if @@ -1647,9 +1710,9 @@ subroutine register_interior_state_variable(self, id, name, units, long_name, & class (type_base_model), intent(inout) :: self type (type_state_variable_id), intent(inout), target :: id character(len=*), intent(in) :: name, long_name, units - real(rk), intent(in), optional :: initial_value,vertical_movement,specific_light_extinction - real(rk), intent(in), optional :: minimum, maximum,missing_value,background_value - logical, intent(in), optional :: no_precipitation_dilution,no_river_dilution + real(rk), intent(in), optional :: initial_value, vertical_movement, specific_light_extinction + real(rk), intent(in), optional :: minimum, maximum, missing_value, background_value + logical, intent(in), optional :: no_precipitation_dilution, no_river_dilution class (type_base_standard_variable), intent(in), optional :: standard_variable integer, intent(in), optional :: presence @@ -1678,7 +1741,6 @@ subroutine register_source(self, link, sms_id, source) source_ = source_do 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-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) @@ -1698,7 +1760,6 @@ subroutine register_surface_flux(self, link, surface_flux_id, source) source_ = source_do_surface 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-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, & @@ -1719,7 +1780,6 @@ subroutine register_bottom_flux(self, link, bottom_flux_id, source) source_ = source_do_bottom 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-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, & @@ -1742,11 +1802,8 @@ subroutine register_movement(self, link, movement_id, vertical_movement) 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-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 - movement_id%link%target%write_operator = operator_add - end if + output=output_none, write_index=movement_id%sum_index, link=movement_id%link, source=source_get_vertical_movement) + movement_id%link%target%write_operator = operator_add link2 => link%target%movement_list%append(movement_id%link%target, movement_id%link%target%name) end subroutine register_movement @@ -1761,7 +1818,6 @@ subroutine register_surface_source(self, link, sms_id, source) source_ = source_do_surface 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-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, & @@ -1782,7 +1838,6 @@ subroutine register_bottom_source(self, link, sms_id, source) source_ = source_do_bottom 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-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, & @@ -1834,7 +1889,7 @@ subroutine add_variable(self, variable, name, units, long_name, missing_value, m initial_value, background_value, fill_value, standard_variable, presence, output, source, & act_as_state_variable, read_index, state_index, write_index, background, link) class (type_base_model), target,intent(inout) :: self - type (type_internal_variable),pointer :: variable + type (type_internal_variable), pointer :: variable character(len=*), target,intent(in) :: name character(len=*), intent(in), optional :: long_name, units real(rk), intent(in), optional :: minimum, maximum, missing_value @@ -2510,54 +2565,25 @@ end subroutine register_named_global_dependency subroutine register_interior_expression_dependency(self, id, expression) class (type_base_model), intent(inout) :: self type (type_dependency_id), target, intent(inout) :: id - class (type_interior_expression), intent(in) :: expression + class (type_interior_expression), pointer :: expression - class (type_interior_expression), allocatable :: copy + class (type_coupling_task), pointer :: base_coupling - allocate(copy, source=expression) - copy%out => id%index - call self%register_dependency(id, copy%output_name, '', copy%output_name) - copy%output_name = id%link%target%name - - call register_expression(self,copy) - deallocate(copy) + call self%register_dependency(id, expression%output_name, '', expression%output_name) + base_coupling => expression + call self%request_coupling(id%link, base_coupling) end subroutine subroutine register_horizontal_expression_dependency(self, id, expression) class (type_base_model), intent(inout) :: self type (type_horizontal_dependency_id), intent(inout), target :: id - class (type_horizontal_expression), intent(in) :: expression - - class (type_horizontal_expression), allocatable :: copy - - allocate(copy, source=expression) - copy%out => id%horizontal_index - call self%register_dependency(id, copy%output_name, '', copy%output_name) - copy%output_name = id%link%target%name + class (type_horizontal_expression), pointer :: expression - call register_expression(self, copy) - deallocate(copy) - end subroutine - - recursive subroutine register_expression(self, expression) - class (type_base_model), intent(inout) :: self - class (type_expression), intent(in) :: expression - - class (type_expression), pointer :: current - - if (.not. associated(self%first_expression)) then - allocate(self%first_expression, source=expression) - current => self%first_expression - else - current => self%first_expression - do while (associated(current%next)) - current => current%next - end do - allocate(current%next, source=expression) - current => current%next - end if + class (type_coupling_task), pointer :: base_coupling - if (associated(self%parent)) call register_expression(self%parent, expression) + call self%register_dependency(id, expression%output_name, '', expression%output_name) + base_coupling => expression + call self%request_coupling(id%link, base_coupling) end subroutine function get_effective_string(value, default) result(value_) @@ -2677,7 +2703,6 @@ function find_object(self, name, recursive, exact) result(object) object => null() link => self%find_link(name, recursive, exact) if (associated(link)) object => link%target - end function find_object recursive function find_link(self, name, recursive, exact) result(link) @@ -2686,9 +2711,9 @@ recursive function find_link(self, name, recursive, exact) result(link) logical, optional, intent(in) :: recursive, exact type (type_link), pointer :: link - integer :: n - logical :: recursive_eff, exact_eff - class (type_base_model),pointer :: current + integer :: n + logical :: recursive_eff, exact_eff + class (type_base_model), pointer :: current link => null() @@ -2785,10 +2810,9 @@ function find_model(self, name, recursive) result(found_model) end do end function find_model - function get_aggregate_variable_access(self, standard_variable, access) result(link) + function get_aggregate_variable_access(self, standard_variable) result(link) class (type_base_model), intent(inout) :: self class (type_domain_specific_standard_variable), target :: standard_variable - integer, optional, intent(in) :: access type (type_link), pointer :: link type (type_aggregate_variable_access), pointer :: aggregate_variable_access @@ -2822,31 +2846,36 @@ function get_aggregate_variable_access(self, standard_variable, access) result(l self%first_aggregate_variable_access => aggregate_variable_access end if link => aggregate_variable_access%link - if (present(access)) link%target%fake_state_variable = link%target%fake_state_variable .or. iand(access, access_add_source) /= 0 end function get_aggregate_variable_access subroutine set_dependency_flag(self, source, flag) class (type_internal_variable), intent(inout) :: self - integer, intent(in) :: source + integer, optional, intent(in) :: source integer, intent(in) :: flag - integer :: n + integer :: source_, n type (type_dependency_flag), allocatable :: prev(:) - n = 1 + if (present(source)) then + source_ = source + else + source_ = -1 + end if if (allocated(self%dependency_flags)) then do n = 1, size(self%dependency_flags) - if (self%dependency_flags(n)%source == source) then + if (self%dependency_flags(n)%source == source_) then ! A flag for this source was already set - overwrite it and return immediately self%dependency_flags(n)%flag = flag return end if end do call move_alloc(self%dependency_flags, prev) + else + n = 1 end if allocate(self%dependency_flags(n)) if (n > 1) self%dependency_flags(:n - 1) = prev - self%dependency_flags(n)%source = source + self%dependency_flags(n)%source = source_ self%dependency_flags(n)%flag = flag end subroutine @@ -2983,18 +3012,21 @@ end subroutine abstract_model_factory_finalize function coupling_task_resolve(self) result(link) class (type_coupling_task), intent(inout) :: self type (type_link), pointer :: link + link => null() end function function link_coupling_task_resolve(self) result(link) class (type_link_coupling_task), intent(inout) :: self type (type_link), pointer :: link + link => self%target_link end function subroutine coupling_task_list_remove(self, task) class (type_coupling_task_list), intent(inout) :: self class (type_coupling_task), pointer :: task + if (associated(task%previous)) then task%previous%next => task%next else @@ -3057,6 +3089,7 @@ end subroutine coupling_task_list_add character(len=32) function source2string(source) integer, intent(in) :: source + select case (source) case (source_unknown); source2string = 'unknown' case (source_state); source2string = 'state' @@ -3077,6 +3110,7 @@ character(len=32) function source2string(source) case (source_get_light_extinction); source2string = 'get_light_extinction' case (source_get_drag); source2string = 'get_drag' case (source_get_albedo); source2string = 'get_albedo' + case (source_global); source2string = 'global' case default write (source2string,'(i0)') source end select @@ -3084,6 +3118,7 @@ end function source2string character(len=32) function domain2string(domain) integer, intent(in) :: domain + select case (domain) case (domain_interior); domain2string = 'interior' case (domain_surface); domain2string = 'surface' @@ -3243,6 +3278,7 @@ subroutine variable_list_finalize(self) function settings_create_child(self) result(child) class (type_fabm_settings), intent(in) :: self class (type_settings), pointer :: child + allocate(type_fabm_settings::child) end function settings_create_child diff --git a/src/fabm_work.F90 b/src/fabm_work.F90 index 1a6e434c..fe504a6c 100644 --- a/src/fabm_work.F90 +++ b/src/fabm_work.F90 @@ -13,80 +13,16 @@ module fabm_work source_check_state, source_check_bottom_state, source_check_surface_state, & domain_interior, domain_surface, domain_bottom, domain_horizontal use fabm_job, only: type_task, type_call + use fabm_global_types use fabm_driver, only: driver implicit none private - public type_domain, type_catalog, type_store, type_cache_fill_values + public type_cache_fill_values public cache_create, cache_pack, cache_unpack - public process_interior_slice, process_horizontal_slice, process_vertical_slice - - ! -------------------------------------------------------------------------- - ! Derived type for model domain - ! (spatial extent, masks, indices of surface and bottom layers) - ! -------------------------------------------------------------------------- - - type type_domain - ! Information about the model domain - integer :: shape(_FABM_DIMENSION_COUNT_) - integer :: start(_FABM_DIMENSION_COUNT_) - integer :: stop(_FABM_DIMENSION_COUNT_) - integer :: horizontal_shape(_HORIZONTAL_DIMENSION_COUNT_) - -#ifdef _HAS_MASK_ -# ifndef _FABM_HORIZONTAL_MASK_ - _FABM_MASK_TYPE_, pointer _ATTRIBUTES_GLOBAL_ :: mask => null() -# endif - _FABM_MASK_TYPE_, pointer _ATTRIBUTES_GLOBAL_HORIZONTAL_ :: mask_hz => null() -#endif - -#ifdef _FABM_DEPTH_DIMENSION_INDEX_ -# if _FABM_BOTTOM_INDEX_==-1 - integer, pointer _ATTRIBUTES_GLOBAL_HORIZONTAL_ :: bottom_indices => null() -# endif -#endif - end type - - ! -------------------------------------------------------------------------- - ! Derived types for catalog with pointers to all available fields - ! -------------------------------------------------------------------------- - - type type_interior_data_pointer - real(rke), pointer _ATTRIBUTES_GLOBAL_ :: p => null() - end type - - type type_horizontal_data_pointer - real(rke), pointer _ATTRIBUTES_GLOBAL_HORIZONTAL_ :: p => null() - end type - - type type_scalar_data_pointer - real(rke), pointer :: p => null() - end type - - type type_catalog - type (type_interior_data_pointer), allocatable :: interior(:) - type (type_horizontal_data_pointer), allocatable :: horizontal(:) - type (type_scalar_data_pointer), allocatable :: scalar(:) - integer, allocatable :: interior_sources(:) - integer, allocatable :: horizontal_sources(:) - integer, allocatable :: scalar_sources(:) - end type - - ! -------------------------------------------------------------------------- - ! Derived type for variable store - ! (spatially explicit model outputs needed by other BGC modules or host) - ! -------------------------------------------------------------------------- - - type type_store - real(rke), allocatable _DIMENSION_GLOBAL_PLUS_1_ :: interior - real(rke), allocatable _DIMENSION_GLOBAL_HORIZONTAL_PLUS_1_ :: horizontal - real(rke), allocatable :: interior_fill_value(:) - real(rke), allocatable :: horizontal_fill_value(:) - real(rke), allocatable :: interior_missing_value(:) - real(rke), allocatable :: horizontal_missing_value(:) - end type + public process_interior_slice, process_horizontal_slice, process_vertical_slice, process_global ! -------------------------------------------------------------------------- ! Derived type for fill/missing values of cache entries @@ -825,6 +761,24 @@ subroutine process_vertical_slice(task, domain, catalog, cache_fill_values, stor end subroutine process_vertical_slice + subroutine process_global(task, catalog _POSTARG_LOCATION_RANGE_, time) + type (type_task), intent(in) :: task + type (type_catalog), intent(in) :: catalog + _DECLARE_ARGUMENTS_LOCATION_RANGE_ + real(rke), optional, intent(in) :: time + + integer :: icall + + do icall = 1, size(task%calls) + if (task%calls(icall)%active) then + select type (model => task%calls(icall)%model) + class is (type_global_model) + call model%update(catalog _POSTARG_LOCATION_RANGE_, time) + end select + end if + end do + end subroutine process_global + subroutine invalidate_interior_call_output(call_node, cache) use fabm_graph, only: type_output_variable_set_node diff --git a/src/models/bb/gas/gas.F90 b/src/models/bb/gas/gas.F90 index 700957e6..d3d5e4a2 100644 --- a/src/models/bb/gas/gas.F90 +++ b/src/models/bb/gas/gas.F90 @@ -28,8 +28,6 @@ subroutine initialize(self, configunit) class (type_bb_gas), intent(inout), target :: self integer, intent(in) :: configunit - call self%register_implemented_routines((/source_do_surface/)) - call self%get_parameter(self%A, 'A', '', 'coefficient A for Schmidt number polynomial (Table 1 in https://doi.org/10.5194/gmd-10-2169-2017)') call self%get_parameter(self%B, 'B', '', 'coefficient B for Schmidt number polynomial (Table 1 in https://doi.org/10.5194/gmd-10-2169-2017)') call self%get_parameter(self%C, 'C', '', 'coefficient C for Schmidt number polynomial (Table 1 in https://doi.org/10.5194/gmd-10-2169-2017)') diff --git a/src/models/examples/mean.F90 b/src/models/examples/mean.F90 index 540d1709..81415b36 100644 --- a/src/models/examples/mean.F90 +++ b/src/models/examples/mean.F90 @@ -47,7 +47,7 @@ subroutine do(self,_ARGUMENTS_DO_) _LOOP_BEGIN_ _GET_(self%id_temp_tempmean, temp) - _SET_DIAGNOSTIC_(self%id_temp_tempmean_diag ,temp) + _SET_DIAGNOSTIC_(self%id_temp_tempmean_diag, temp) _LOOP_END_ end subroutine do diff --git a/src/models/su/mixo.F90 b/src/models/su/mixo.F90 index 889a92a2..f942f9fc 100644 --- a/src/models/su/mixo.F90 +++ b/src/models/su/mixo.F90 @@ -400,6 +400,8 @@ subroutine initialize(self,configunit) call self%register_dependency(self%id_mPStotdep,'mPStot', 'gC (gC)-1 (day)-1', 'enables reading the previous value of the diagnostic mPStot') call self%register_dependency(self%id_avgCu, temporal_mean(self%id_mCudep, period=86400._rk,resolution=3600._rk,missing_value=0.0_rk)) call self%register_dependency(self%id_avgPStot, temporal_mean(self%id_mPStotdep, period=86400._rk,resolution=3600._rk,missing_value=0.0_rk)) + call set_dependency_flag(self%id_avgCu%link%target, flag=dependency_flag_stale) + call set_dependency_flag(self%id_avgPStot%link%target, flag=dependency_flag_stale) call self%register_dependency(self%id_PAR, standard_variables%downwelling_photosynthetic_radiative_flux) ! local PAR call self%register_dependency(self%id_ETW, standard_variables%temperature) diff --git a/src/pyfabm/__init__.py b/src/pyfabm/__init__.py index 2cd61886..b3e1440a 100644 --- a/src/pyfabm/__init__.py +++ b/src/pyfabm/__init__.py @@ -165,20 +165,8 @@ def get_lib(name: str) -> FABMDLL: # Access to model objects (variables, parameters, dependencies, couplings, # model instances) - lib.get_counts.argtypes = [ - ctypes.c_void_p, - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - ] + ncat = 12 if hasattr(lib, "get_variable_part_of_state") else 11 + lib.get_counts.argtypes = [ctypes.c_void_p] + [ctypes.POINTER(ctypes.c_int)] * ncat lib.get_counts.restype = None lib.get_variable_metadata.argtypes = [ ctypes.c_void_p, @@ -198,6 +186,15 @@ def get_lib(name: str) -> FABMDLL: ctypes.c_int, ] lib.set_variable_save.restype = None + try: + lib.get_variable_part_of_state.argtypes = [ + ctypes.c_void_p, + ctypes.c_int, + ctypes.c_int, + ] + lib.get_variable_part_of_state.restype = ctypes.c_int + except AttributeError: + setattr(lib, "get_variable_part_of_state", lambda p, t, i: 0) lib.get_variable.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int] lib.get_variable.restype = ctypes.c_void_p lib.get_parameter_metadata.argtypes = [ @@ -391,6 +388,11 @@ def get_lib(name: str) -> FABMDLL: lib.get_interior_diagnostic_data.restype = ctypes.POINTER(lib.dtype) lib.get_horizontal_diagnostic_data.argtypes = [ctypes.c_void_p, ctypes.c_int] lib.get_horizontal_diagnostic_data.restype = ctypes.POINTER(lib.dtype) + try: + lib.get_scalar_diagnostic_data.argtypes = [ctypes.c_void_p, ctypes.c_int] + lib.get_scalar_diagnostic_data.restype = ctypes.POINTER(lib.dtype) + except AttributeError: + pass lib.require_data.argtypes = [ctypes.c_void_p, ctypes.c_void_p] lib.require_data.restype = None lib.get_standard_variable_data.argtypes = [ @@ -427,7 +429,7 @@ def get_lib(name: str) -> FABMDLL: lib.check_state.restype = ctypes.c_int # Routine for getting git repository version information. - lib.get_version.argtypes = (ctypes.c_int, ctypes.c_char_p) + lib.get_version.argtypes = [ctypes.c_int, ctypes.c_char_p] lib.get_version.restype = None lib.save_settings.argtypes = [ctypes.c_void_p, ctypes.c_char_p, ctypes.c_int] @@ -478,6 +480,7 @@ def log(msg: str): INTERIOR_DEPENDENCY = 7 HORIZONTAL_DEPENDENCY = 8 SCALAR_DEPENDENCY = 9 +SCALAR_DIAGNOSTIC_VARIABLE = 10 class DataType(enum.IntEnum): @@ -627,9 +630,16 @@ def __init__( path: Optional[str] = None, ): self.model = model + + #: Short variable name (alphanumeric characters and underscores only) self.name = name + + #: Units self.units = units + self.units_unicode = None if units is None else createPrettyUnit(units) + + #: Long name self.long_name = long_name or name self.path = path or name @@ -645,7 +655,7 @@ def output_name(self) -> str: @property def options(self) -> Optional[Sequence]: """Collection of values that this variable can take. - `None` if the variable is not limited to any particular value.""" + ``None`` if the variable is not limited to any particular value.""" return None def __repr__(self) -> str: @@ -779,11 +789,11 @@ def __init__( model: "Model", variable_pointer: ctypes.c_void_p, index: int, - horizontal: bool, + type: int, ): super().__init__(model, variable_pointer) self._data: Optional[np.ndarray] = None - self._horizontal = horizontal + self._type = type self._index = index + 1 @property @@ -804,15 +814,20 @@ def save(self) -> bool: @save.setter def save(self, value: bool): - vartype = ( - HORIZONTAL_DIAGNOSTIC_VARIABLE - if self._horizontal - else INTERIOR_DIAGNOSTIC_VARIABLE - ) self.model.fabm.set_variable_save( - self.model.pmodel, vartype, self._index, 1 if value else 0 + self.model.pmodel, self._type, self._index, 1 if value else 0 ) + @property + def part_of_state(self) -> bool: + """Whether the this diagnostic is part of the model state, + (that is, its current value impacts future source terms) + and therefore needs to be included in restarts""" + value: int = self.model.fabm.get_variable_part_of_state( + self.model.pmodel, self._type, self._index + ) + return value != 0 + class Parameter(Variable): def __init__( @@ -1117,6 +1132,7 @@ def __init__( self.bottom_state_variables: NamedObjectList[StateVariable] = NamedObjectList() self.interior_diagnostic_variables: NamedObjectList[DiagnosticVariable] = NamedObjectList() self.horizontal_diagnostic_variables: NamedObjectList[DiagnosticVariable] = NamedObjectList() + self.scalar_diagnostic_variables: NamedObjectList[DiagnosticVariable] = NamedObjectList() self.conserved_quantities: NamedObjectList[ConservedQuantity] = NamedObjectList() self.parameters: NamedObjectList[Parameter] = NamedObjectList() self.interior_dependencies: NamedObjectList[Dependency] = NamedObjectList() @@ -1304,26 +1320,33 @@ def _update_configuration( nstate_bottom = ctypes.c_int() ndiag_interior = ctypes.c_int() ndiag_horizontal = ctypes.c_int() + ndiag_scalar = ctypes.c_int(0) ndependencies_interior = ctypes.c_int() ndependencies_horizontal = ctypes.c_int() ndependencies_scalar = ctypes.c_int() nconserved = ctypes.c_int() nparameters = ctypes.c_int() ncouplings = ctypes.c_int() - self.fabm.get_counts( - self.pmodel, + args = [] + args += [ ctypes.byref(nstate_interior), ctypes.byref(nstate_surface), ctypes.byref(nstate_bottom), - ctypes.byref(ndiag_interior), - ctypes.byref(ndiag_horizontal), + ] + args += [ctypes.byref(ndiag_interior), ctypes.byref(ndiag_horizontal)] + if len(self.fabm.get_counts.argtypes) == 13: + args += [ctypes.byref(ndiag_scalar)] + args += [ ctypes.byref(ndependencies_interior), ctypes.byref(ndependencies_horizontal), ctypes.byref(ndependencies_scalar), + ] + args += [ ctypes.byref(nconserved), ctypes.byref(nparameters), ctypes.byref(ncouplings), - ) + ] + self.fabm.get_counts(self.pmodel, *args) # Allocate memory for state variable values, and send ctypes.pointer to # this memory to FABM. @@ -1373,6 +1396,7 @@ def _update_configuration( self.bottom_state_variables.clear() self.interior_diagnostic_variables.clear() self.horizontal_diagnostic_variables.clear() + self.scalar_diagnostic_variables.clear() self.conserved_quantities.clear() self.parameters.clear() self.interior_dependencies.clear() @@ -1398,14 +1422,19 @@ def _update_configuration( self.pmodel, INTERIOR_DIAGNOSTIC_VARIABLE, i + 1 ) self.interior_diagnostic_variables._data.append( - DiagnosticVariable(self, ptr, i, False) + DiagnosticVariable(self, ptr, i, INTERIOR_DIAGNOSTIC_VARIABLE) ) for i in range(ndiag_horizontal.value): ptr = self.fabm.get_variable( self.pmodel, HORIZONTAL_DIAGNOSTIC_VARIABLE, i + 1 ) self.horizontal_diagnostic_variables._data.append( - DiagnosticVariable(self, ptr, i, True) + DiagnosticVariable(self, ptr, i, HORIZONTAL_DIAGNOSTIC_VARIABLE) + ) + for i in range(ndiag_scalar.value): + ptr = self.fabm.get_variable(self.pmodel, SCALAR_DIAGNOSTIC_VARIABLE, i + 1) + self.scalar_diagnostic_variables._data.append( + DiagnosticVariable(self, ptr, i, SCALAR_DIAGNOSTIC_VARIABLE) ) for i in range(ndependencies_interior.value): ptr = self.fabm.get_variable(self.pmodel, INTERIOR_DEPENDENCY, i + 1) @@ -1483,7 +1512,9 @@ def _update_configuration( + self.bottom_state_variables ) self.diagnostic_variables = ( - self.interior_diagnostic_variables + self.horizontal_diagnostic_variables + self.interior_diagnostic_variables + + self.horizontal_diagnostic_variables + + self.scalar_diagnostic_variables ) self.dependencies = ( self.interior_dependencies @@ -1705,6 +1736,13 @@ def process_dependencies(dependencies: Sequence[Dependency]): variable._data = arr.view(dtype=self.fabm.numpy_dtype) else: variable._data = None + for i, variable in enumerate(self.scalar_diagnostic_variables): + pdata = self.fabm.get_scalar_diagnostic_data(self.pmodel, i + 1) + if pdata: + arr = np.ctypeslib.as_array(pdata, ()) + variable._data = arr.view(dtype=self.fabm.numpy_dtype) + else: + variable._data = None return ready checkReady = start diff --git a/src/pyfabm/gui_qt.py b/src/pyfabm/gui_qt.py index 47bb653d..bb59aacf 100644 --- a/src/pyfabm/gui_qt.py +++ b/src/pyfabm/gui_qt.py @@ -344,7 +344,7 @@ def setData( QtCore.Qt.ItemDataRole.CheckStateRole, ): assert isinstance(value, (float, int, bool, str)) - entry = cast(Entry, index.internalPointer()) + entry: Entry = index.internalPointer() data = entry.object assert isinstance(data, pyfabm.Variable) data.value = value diff --git a/testcases/fabm-builtin.yaml b/testcases/fabm-builtin.yaml index 86230336..d109d199 100644 --- a/testcases/fabm-builtin.yaml +++ b/testcases/fabm-builtin.yaml @@ -123,4 +123,31 @@ instances: maximum_depth: 20.0 coupling: source: tracer/c - \ No newline at end of file + interior_temporal_mean: + model: interior_temporal_mean + parameters: + window: 86400 + n: 24 + coupling: + source: tracer/c + surface_temporal_mean: + model: surface_temporal_mean + parameters: + window: 86400 + n: 24 + coupling: + source: surface_layer/result + bottom_temporal_mean: + model: bottom_temporal_mean + parameters: + window: 86400 + n: 24 + coupling: + source: bottom_layer/result + surface_temporal_maximum: + model: surface_temporal_maximum + parameters: + window: 86400 + n: 24 + coupling: + source: surface_layer/result diff --git a/util/algorithms/running_max.ipynb b/util/algorithms/running_max.ipynb index a670aac1..662e7418 100644 --- a/util/algorithms/running_max.ipynb +++ b/util/algorithms/running_max.ipynb @@ -2,223 +2,248 @@ "cells": [ { "cell_type": "markdown", + "metadata": {}, "source": [ - "# Prototype: an efficient, approximate running maximum filter\r\n", - "\r\n", - "The window size (units of time) is denoted as $p$.\r\n", - "\r\n", - "Constraints:\r\n", - "* The filter uses a prescribed number of bins, $n$, to store the history.\r\n", - " Typically $n$ << the number of time steps the model would take to cover the window.\r\n", - " For instance, a 1-year running mean may use just 1 bin. In that case, the maximum during very first year will be set to \"missing value\".\r\n", - " In the entire subsequent year, the maximum will be equal to the maximum over the first year. And so on.\r\n", - " By definition, this means the running maximum will be an approximation that ignores all recent history (of duration window size / # bins).\r\n", + "# Prototype: an efficient, approximate running maximum filter\n", + "\n", + "The window size (units of time) is denoted as $p$.\n", + "\n", + "Constraints:\n", + "* The filter uses a prescribed number of bins, $n$, to store the history.\n", + " Typically $n$ << the number of time steps the model would take to cover the window.\n", + " For instance, a 1-year running mean may use just 1 bin. In that case, the maximum during very first year will be set to \"missing value\".\n", + " In the entire subsequent year, the maximum will be equal to the maximum over the first year. And so on.\n", + " By definition, this means the running maximum will be an approximation that ignores all recent history (of duration window size / # bins).\n", "* We have no advance knowledge about the model time step $\\Delta t$, which may be variable." - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, - "source": [ - "import numpy\r\n", - "from matplotlib import pyplot\r\n", - "\r\n", - "p = 1. # window size time units\r\n", - "delta_t = 1. / 24. / 6. / 365 # model time step in time units\r\n", - "duration = 5. # simulation duration in time units\r\n", - "n = 1 # number of bins for history (covering period p)\r\n", - "missing_value = -2. # value to return while the simulation has not covered 1 window size yet" - ], + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "import numpy as np\n", + "from matplotlib import pyplot\n", + "\n", + "p = 1.0 # window size time units\n", + "delta_t = 1.0 / 24.0 / 6.0 / 365 # model time step in time units\n", + "duration = 5.0 # simulation duration in time units\n", + "n = 1 # number of bins for history (covering period p)\n", + "missing_value = -2.0 # value to return while the simulation has not covered 1 window size yet" + ] }, { "cell_type": "code", "execution_count": null, - "source": [ - "class Filter:\r\n", - " def __init__(self):\r\n", - " self.history = numpy.zeros((n,))\r\n", - " self.previous_time = 0.\r\n", - " self.previous_value = 0\r\n", - " self.ibin = -1\r\n", - " self.bin_end_time = 0.\r\n", - " self.max = missing_value\r\n", - " self.complete = False\r\n", - "\r\n", - " def __call__(self, now: float, value: float) -> float:\r\n", - " if self.ibin == -1:\r\n", - " self.bin_end_time = now + p / n\r\n", - " self.ibin = 0\r\n", - "\r\n", - " while now >= self.bin_end_time:\r\n", - " # Interpolate to value at right bin time\r\n", - " w = (self.bin_end_time - self.previous_time) / (now - self.previous_time)\r\n", - " self.previous_time = self.bin_end_time\r\n", - " self.previous_value += w * (value - self.previous_value)\r\n", - "\r\n", - " # Increment the bin we are completing (history[ibin]) and mean\r\n", - " #bin_end_value = (1 - w) * self.previous_value + w * value\r\n", - " #self.history[self.ibin] += 0.5 * dt * (self.previous_value + bin_end_value) / p\r\n", - " self.history[self.ibin] = numpy.maximum(self.history[self.ibin], self.previous_value)\r\n", - " self.complete = self.complete or self.ibin == n - 1\r\n", - " if self.complete:\r\n", - " self.max = numpy.max(self.history[:, ...], axis=0)\r\n", - " self.ibin = 0 if self.ibin == n - 1 else self.ibin + 1\r\n", - " self.history[self.ibin] = self.previous_value\r\n", - "\r\n", - " self.bin_end_time += p / n\r\n", - "\r\n", - " # Update the maximum of the current bin\r\n", - " self.history[self.ibin] = numpy.maximum(self.history[self.ibin], value)\r\n", - "\r\n", - " # Store current time and value to enable linear interpolation in subsequent call.\r\n", - " self.previous_time = now\r\n", - " self.previous_value = value\r\n", - " return self.max" - ], + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "class Filter:\n", + " def __init__(self):\n", + " self.history = np.zeros((n,))\n", + " self.previous_time = 0.0\n", + " self.previous_value = 0\n", + " self.ibin = -1\n", + " self.bin_end_time = 0.0\n", + " self.max = missing_value\n", + " self.complete = False\n", + "\n", + " def __call__(self, now: float, value: float) -> float:\n", + " if self.ibin == -1:\n", + " self.bin_end_time = now + p / n\n", + " self.ibin = 0\n", + "\n", + " while now >= self.bin_end_time:\n", + " # Interpolate to value at right bin time\n", + " w = (self.bin_end_time - self.previous_time) / (now - self.previous_time)\n", + " self.previous_time = self.bin_end_time\n", + " self.previous_value += w * (value - self.previous_value)\n", + "\n", + " # Increment the bin we are completing (history[ibin]) and mean\n", + " # bin_end_value = (1 - w) * self.previous_value + w * value\n", + " # self.history[self.ibin] += 0.5 * dt * (self.previous_value + bin_end_value) / p\n", + " self.history[self.ibin] = np.maximum(\n", + " self.history[self.ibin], self.previous_value\n", + " )\n", + " self.complete = self.complete or self.ibin == n - 1\n", + " if self.complete:\n", + " self.max = np.max(self.history[:, ...], axis=0)\n", + " self.ibin = 0 if self.ibin == n - 1 else self.ibin + 1\n", + " self.history[self.ibin] = self.previous_value\n", + "\n", + " self.bin_end_time += p / n\n", + "\n", + " # Update the maximum of the current bin\n", + " self.history[self.ibin] = np.maximum(self.history[self.ibin], value)\n", + "\n", + " # Store current time and value to enable linear interpolation in subsequent call.\n", + " self.previous_time = now\n", + " self.previous_value = value\n", + " return self.max" + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ - "# Calculate and plot variable for which to compute the running maximum\r\n", - "times = numpy.arange(0, duration, delta_t)\r\n", - "values = numpy.sin(2 * numpy.pi * times / 3)\r\n", - "\r\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\r\n", - "ax.plot(times, values)\r\n", + "# Calculate and plot variable for which to compute the running maximum\n", + "times = np.arange(0, duration, delta_t)\n", + "values = np.sin(2 * np.pi * times / 3)\n", + "\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "ax.plot(times, values)\n", "ax.grid()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ - "# Compute and plot the running maximum\r\n", - "filter = Filter()\r\n", - "filtered = numpy.empty_like(values)\r\n", - "for i, (time, value) in enumerate(zip(times, values)):\r\n", - " filtered[i] = filter(time, value)\r\n", - "\r\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\r\n", - "ax.plot(times, values)\r\n", - "ax.plot(times, filtered)\r\n", + "# Compute and plot the running maximum\n", + "filter = Filter()\n", + "filtered = np.empty_like(values)\n", + "for i, (time, value) in enumerate(zip(times, values)):\n", + " filtered[i] = filter(time, value)\n", + "\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "ax.plot(times, values)\n", + "ax.plot(times, filtered)\n", "ax.grid()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, - "source": [ - "# Compare running maximum with analytical solution\r\n", - "# (this requires the window size to be a multiple of the model time step)\r\n", - "assert abs(p % delta_t) < 1e-15, 'Window size %s is not a multiple of the model time step %s. Residual: %s' % (p, delta_t, p % delta_t)\r\n", - "nstep = int(round(p / delta_t))\r\n", - "analytical = numpy.full_like(values, missing_value)\r\n", - "for i, (time, value) in enumerate(zip(times, values)):\r\n", - " if i >= nstep:\r\n", - " analytical[i] = values[i - nstep:i + 1].max()\r\n", - "\r\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\r\n", - "ax.plot(times, filtered - analytical)\r\n", - "ax.grid()\r\n" - ], + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "# Compare running maximum with analytical solution\n", + "# (this requires the window size to be a multiple of the model time step)\n", + "assert (\n", + " abs(p % delta_t) < 1e-15\n", + "), \"Window size %s is not a multiple of the model time step %s. Residual: %s\" % (\n", + " p,\n", + " delta_t,\n", + " p % delta_t,\n", + ")\n", + "nstep = int(round(p / delta_t))\n", + "analytical = np.full_like(values, missing_value)\n", + "for i, (time, value) in enumerate(zip(times, values)):\n", + " if i >= nstep:\n", + " analytical[i] = values[i - nstep : i + 1].max()\n", + "\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "ax.plot(times, filtered - analytical)\n", + "ax.grid()" + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ - "# Now try with randomly varying time step\r\n", - "filter = Filter()\r\n", - "time = 0.\r\n", - "rtimes, rvalues, rfiltered = [], [], []\r\n", - "while time < duration:\r\n", - " value = numpy.sin(2 * numpy.pi * time / 3)\r\n", - " rtimes.append(time)\r\n", - " rvalues.append(value)\r\n", - " rfiltered.append(filter(time, value))\r\n", - " dt = 2 * delta_t * numpy.random.rand()\r\n", - " time += dt\r\n", - "\r\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\r\n", - "ax.plot(rtimes, rvalues)\r\n", - "ax.plot(rtimes, rfiltered)\r\n", + "# Now try with randomly varying time step\n", + "filter = Filter()\n", + "time = 0.0\n", + "rtimes, rvalues, rfiltered = [], [], []\n", + "while time < duration:\n", + " value = np.sin(2 * np.pi * time / 3)\n", + " rtimes.append(time)\n", + " rvalues.append(value)\n", + " rfiltered.append(filter(time, value))\n", + " dt = 2 * delta_t * np.random.rand()\n", + " time += dt\n", + "\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "ax.plot(rtimes, rvalues)\n", + "ax.plot(rtimes, rfiltered)\n", "ax.grid()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ - "# Now try the running maximum filter in FABM itself\r\n", - "import pyfabm\r\n", - "config = {'instances': {'max': {'model': 'surface_temporal_maximum', 'parameters': {'window': p, 'n': n, 'missing_value': missing_value}}}}\r\n", - "m = pyfabm.Model(config)\r\n", - "invar = m.findDependency('max/source')\r\n", - "outvar = m.findDiagnosticVariable('max/max')\r\n", - "invar.value = missing_value\r\n", - "m.cell_thickness = 1\r\n", - "m.start()\r\n", - "fabm_fitered = numpy.empty_like(values)\r\n", - "for i, (time, value) in enumerate(zip(times, values)):\r\n", - " invar.value = value\r\n", - " m.getRates(time, surface=True)\r\n", - " fabm_fitered[i] = outvar.value\r\n", - "\r\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\r\n", - "ax.plot(times, values)\r\n", - "ax.plot(times, fabm_fitered)\r\n", - "ax.grid()\r\n", - "\r\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\r\n", - "ax.plot(times, fabm_fitered - analytical)\r\n", - "ax.grid()\r\n" - ], + "# Now try the running maximum filter in FABM itself\n", + "import pyfabm\n", + "\n", + "config = {\n", + " \"instances\": {\n", + " \"max\": {\n", + " \"model\": \"surface_temporal_maximum\",\n", + " \"parameters\": {\"window\": p, \"n\": n, \"missing_value\": missing_value},\n", + " }\n", + " }\n", + "}\n", + "m = pyfabm.Model(config)\n", + "invar = m.dependencies[\"max/source\"]\n", + "outvar = m.diagnostic_variables[\"max/max\"]\n", + "invar.value = missing_value\n", + "m.cell_thickness = 1.0\n", + "m.start()\n", + "fabm_filtered = np.empty_like(values)\n", + "for i, (time, value) in enumerate(zip(times, values)):\n", + " invar.value = value\n", + " m.get_sources(time)\n", + " fabm_filtered[i] = outvar.value\n", + "\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "ax.plot(times, values)\n", + "ax.plot(times, fabm_filtered)\n", + "ax.grid()\n", + "\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "ax.plot(times, fabm_filtered - analytical)\n", + "ax.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "fig, ax = pyplot.subplots()\n", + "ax.plot(filtered, fabm_filtered, \".\")\n", + "ax.grid()" + ] }, { "cell_type": "code", "execution_count": null, - "source": [], + "metadata": {}, "outputs": [], - "metadata": {} + "source": [] } ], "metadata": { - "orig_nbformat": 4, + "kernelspec": { + "display_name": "base", + "language": "python", + "name": "python3" + }, "language_info": { - "name": "python", - "version": "3.8.8", - "mimetype": "text/x-python", "codemirror_mode": { "name": "ipython", "version": 3 }, - "pygments_lexer": "ipython3", + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", "nbconvert_exporter": "python", - "file_extension": ".py" - }, - "kernelspec": { - "name": "python3", - "display_name": "Python 3.8.8 64-bit ('base': conda)" + "pygments_lexer": "ipython3", + "version": "3.9.21" }, - "interpreter": { - "hash": "97ae724bfa85b9b34df7982b8bb8c7216f435b92902d749e4263f71162bea840" - } + "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 -} \ No newline at end of file +} diff --git a/util/algorithms/running_mean.ipynb b/util/algorithms/running_mean.ipynb index 827dae57..b6d6ff6f 100644 --- a/util/algorithms/running_mean.ipynb +++ b/util/algorithms/running_mean.ipynb @@ -23,14 +23,14 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy\n", + "import numpy as np\n", "from matplotlib import pyplot\n", "\n", - "p = 1.0 # window size in days\n", - "delta_t = 1.0 / 24.0 / 6.0 # model time step in days\n", - "duration = 30.0 # simulation duration in days\n", - "n = 24 # number of bins for history (covering period p)\n", - "missing_value = -2.0 # value to return while the simulation has not covered 1 window size yet" + "p = 1.0 # window size in days\n", + "delta_t = 1.0 / 24.0 / 6.0 # model time step in days\n", + "duration = 30.0 # simulation duration in days\n", + "n = 24 # number of bins for history (covering period p)\n", + "missing_value = -2.0 # value to return while the simulation has not covered 1 window size yet" ] }, { @@ -40,56 +40,71 @@ "outputs": [], "source": [ "class Filter:\n", - " def __init__(self):\n", - " self.history = numpy.zeros((n + 1,)) # one for each completed bin, plus one more for the bin currently being filled\n", - " self.previous_time = 0.\n", - " self.previous_value = 0\n", - " self.ibin = -1\n", - " self.bin_end_time = 0.\n", - " self.mean_value = missing_value\n", - " self.last_mean = 0.\n", - " self.complete = False\n", - "\n", - " def __call__(self, now: float, value: float) -> float:\n", - " binwidth = p / n\n", - " if self.ibin == -1:\n", - " self.previous_time = now\n", - " self.bin_end_time = now + binwidth\n", - " self.ibin = 0\n", - "\n", - " while now >= self.bin_end_time:\n", - " dt = self.bin_end_time - self.previous_time\n", - "\n", - " # Interpolate to value at right bin time\n", - " w = dt / (now - self.previous_time)\n", - "\n", - " # Increment the bin we are completing (history[ibin]) and mean\n", - " #bin_end_value = (1 - w) * self.previous_value + w * value\n", - " #self.history[self.ibin] += 0.5 * dt * (self.previous_value + bin_end_value) / p\n", - " self.history[self.ibin] += (self.previous_value + 0.5 * w * (value - self.previous_value)) * dt / p\n", - " if self.complete:\n", - " # We already had a complete history (bins covering the full window size). Add the newly full bin, subtract the oldest bin\n", - " self.last_mean += self.history[self.ibin] - self.history[self.ibin + 1 if self.ibin != n else 0]\n", - " elif self.ibin == n - 1:\n", - " # We just completed our history. create the mean by summing all filled bins.\n", - " self.last_mean = self.history[:n, ...].sum(axis=0)\n", - " self.complete = True\n", - " self.ibin = 0 if self.ibin == n else self.ibin + 1\n", - " self.history[self.ibin] = 0.\n", - "\n", - " self.previous_time = self.bin_end_time\n", - " self.previous_value += w * (value - self.previous_value)\n", - " self.bin_end_time += binwidth\n", - "\n", - " # increment current bin (history[ibin])\n", - " self.history[self.ibin] += 0.5 * (self.previous_value + value) / p * (now - self.previous_time)\n", - " if self.complete:\n", - " # we have a complete history - update the mean\n", - " self.mean_value = self.last_mean + self.history[self.ibin] - self.history[self.ibin + 1 if self.ibin != n else 0] * (now - self.bin_end_time + binwidth) / binwidth\n", - "\n", - " self.previous_time = now\n", - " self.previous_value = value\n", - " return self.mean_value" + " def __init__(self):\n", + " self.history = np.zeros(\n", + " (n + 1,)\n", + " ) # one for each completed bin, plus one more for the bin currently being filled\n", + " self.previous_time = 0.0\n", + " self.previous_value = 0\n", + " self.ibin = -1\n", + " self.bin_end_time = 0.0\n", + " self.mean_value = missing_value\n", + " self.last_mean = 0.0\n", + " self.complete = False\n", + "\n", + " def __call__(self, now: float, value: float) -> float:\n", + " binwidth = p / n\n", + " if self.ibin == -1:\n", + " self.previous_time = now\n", + " self.bin_end_time = now + binwidth\n", + " self.ibin = 0\n", + "\n", + " while now >= self.bin_end_time:\n", + " dt = self.bin_end_time - self.previous_time\n", + "\n", + " # Interpolate to value at right bin time\n", + " w = dt / (now - self.previous_time)\n", + "\n", + " # Increment the bin we are completing (history[ibin]) and mean\n", + " # bin_end_value = (1 - w) * self.previous_value + w * value\n", + " # self.history[self.ibin] += 0.5 * dt * (self.previous_value + bin_end_value) / p\n", + " self.history[self.ibin] += (\n", + " (self.previous_value + 0.5 * w * (value - self.previous_value)) * dt / p\n", + " )\n", + " if self.complete:\n", + " # We already had a complete history (bins covering the full window size). Add the newly full bin, subtract the oldest bin\n", + " self.last_mean += (\n", + " self.history[self.ibin]\n", + " - self.history[self.ibin + 1 if self.ibin != n else 0]\n", + " )\n", + " elif self.ibin == n - 1:\n", + " # We just completed our history. create the mean by summing all filled bins.\n", + " self.last_mean = self.history[:n, ...].sum(axis=0)\n", + " self.complete = True\n", + " self.ibin = 0 if self.ibin == n else self.ibin + 1\n", + " self.history[self.ibin] = 0.0\n", + "\n", + " self.previous_time = self.bin_end_time\n", + " self.previous_value += w * (value - self.previous_value)\n", + " self.bin_end_time += binwidth\n", + "\n", + " # increment current bin (history[ibin])\n", + " self.history[self.ibin] += (\n", + " 0.5 * (self.previous_value + value) / p * (now - self.previous_time)\n", + " )\n", + " if self.complete:\n", + " # we have a complete history - update the mean\n", + " self.mean_value = (\n", + " self.last_mean\n", + " + self.history[self.ibin]\n", + " - self.history[self.ibin + 1 if self.ibin != n else 0]\n", + " * (now - self.bin_end_time + binwidth)\n", + " / binwidth\n", + " )\n", + "\n", + " self.previous_time = now\n", + " self.previous_value = value\n", + " return self.mean_value" ] }, { @@ -99,10 +114,10 @@ "outputs": [], "source": [ "# Calculate and plot variable for which to compute the running mean\n", - "times = numpy.arange(0, duration, delta_t)\n", - "values = numpy.sin(2 * numpy.pi * times)\n", + "times = np.arange(0, duration, delta_t)\n", + "values = np.sin(2 * np.pi * times)\n", "\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", "ax.plot(times, values)\n", "ax.grid()" ] @@ -115,11 +130,11 @@ "source": [ "# Compute and plot the running mean\n", "filter = Filter()\n", - "filtered = numpy.empty_like(values)\n", + "filtered = np.empty_like(values)\n", "for i, (time, value) in enumerate(zip(times, values)):\n", - " filtered[i] = filter(time, value)\n", + " filtered[i] = filter(time, value)\n", "\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", "ax.plot(times, values)\n", "ax.plot(times, filtered)\n", "ax.grid()" @@ -133,17 +148,23 @@ "source": [ "# Compare running mean with analytical solution\n", "# (this requires the window size to be a multiple of the model time step)\n", - "assert abs(p % delta_t) < 1e-15, 'Window size %s is not a multiple of the model time step %s. Residual: %s' % (p, delta_t, p % delta_t)\n", + "assert (\n", + " abs(p % delta_t) < 1e-15\n", + "), \"Window size %s is not a multiple of the model time step %s. Residual: %s\" % (\n", + " p,\n", + " delta_t,\n", + " p % delta_t,\n", + ")\n", "nstep = int(round(p / delta_t))\n", - "analytical = numpy.full_like(values, missing_value)\n", + "analytical = np.full_like(values, missing_value)\n", "for i, (time, value) in enumerate(zip(times, values)):\n", - " if i >= nstep:\n", - " centers = 0.5 * (values[i - nstep:i] + values[i - nstep + 1:i + 1])\n", - " analytical[i] = centers.mean()\n", + " if i >= nstep:\n", + " centers = 0.5 * (values[i - nstep : i] + values[i - nstep + 1 : i + 1])\n", + " analytical[i] = centers.mean()\n", "\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", "ax.plot(times, filtered - analytical)\n", - "ax.grid()\n" + "ax.grid()" ] }, { @@ -154,17 +175,17 @@ "source": [ "# Now try with randomly varying time step\n", "filter = Filter()\n", - "time = 0.\n", + "time = 0.0\n", "rtimes, rvalues, rfiltered = [], [], []\n", "while time < duration:\n", - " value = numpy.sin(2 * numpy.pi * time)\n", - " rtimes.append(time)\n", - " rvalues.append(value)\n", - " rfiltered.append(filter(time, value))\n", - " dt = 2 * delta_t * numpy.random.rand()\n", - " time += dt\n", - "\n", - "fig, ax = pyplot.subplots(figsize=(15, 4))\n", + " value = np.sin(2 * np.pi * time)\n", + " rtimes.append(time)\n", + " rvalues.append(value)\n", + " rfiltered.append(filter(time, value))\n", + " dt = 2 * delta_t * np.random.rand()\n", + " time += dt\n", + "\n", + "fig, ax = pyplot.subplots(figsize=(15, 4))\n", "ax.plot(rtimes, rvalues)\n", "ax.plot(rtimes, rfiltered)\n", "ax.grid()" @@ -188,7 +209,7 @@ " mean_sf=dict(\n", " model=\"surface_temporal_mean\",\n", " parameters=dict(window=p, n=n, missing_value=missing_value),\n", - " )\n", + " ),\n", " )\n", ")\n", "m = pyfabm.Model(config)\n", @@ -200,17 +221,17 @@ "invar2.value = missing_value\n", "m.cell_thickness = 1.0\n", "m.start()\n", - "fabm_fitered = numpy.empty_like(values)\n", + "fabm_filtered = np.empty_like(values)\n", "for i, (time, value) in enumerate(zip(times, values)):\n", " invar.value = value\n", " invar2.value = value\n", - " m.getRates(time, surface=True)\n", + " m.get_sources(time)\n", " assert outvar.value == outvar2.value\n", - " fabm_fitered[i] = outvar.value\n", + " fabm_filtered[i] = outvar.value\n", "\n", "fig, ax = pyplot.subplots(figsize=(15, 4))\n", - "ax.plot(times, fabm_fitered - analytical)\n", - "ax.grid()\n" + "ax.plot(times, fabm_filtered - analytical)\n", + "ax.grid()" ] }, { @@ -222,11 +243,9 @@ } ], "metadata": { - "interpreter": { - "hash": "97ae724bfa85b9b34df7982b8bb8c7216f435b92902d749e4263f71162bea840" - }, "kernelspec": { - "display_name": "Python 3.8.8 64-bit ('base': conda)", + "display_name": "base", + "language": "python", "name": "python3" }, "language_info": { @@ -239,7 +258,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.16" + "version": "3.9.21" }, "orig_nbformat": 4 },