Commit d5c7fbe9 authored by Jorn Bruggeman's avatar Jorn Bruggeman
Browse files

python driver: do not ignore bottom fluxes for pelagic variables; divide...

python driver: do not ignore bottom fluxes for pelagic variables; divide bottom/surface fluxes by column height
parent 15c9a5bb
......@@ -24,13 +24,16 @@
contains
subroutine get_environment_metadata(model,environment_names,environment_units)
subroutine get_environment_metadata(model,environment_names,environment_units,index_column_depth)
type (type_model), intent(inout) :: model
character(len=1024),dimension(:),allocatable,intent(out) :: environment_names,environment_units
integer, intent(out) :: index_column_depth
integer :: n
type (type_link), pointer :: link
index_column_depth = -1
! Get number of environmental dependencies (light, temperature, etc.)
n = 0
link => model%links_postcoupling%first
......@@ -39,8 +42,14 @@
select case (link%target%domain)
case (domain_interior)
if (.not.associated(model%data(link%target%read_indices%pointers(1)%p)%p)) n = n+1
if (index_column_depth==-1 .and. associated(link%target%standard_variables%first)) then
if (link%target%standard_variables%first%p%compare(standard_variables%cell_thickness)) index_column_depth = n
end if
case (domain_bottom,domain_surface,domain_horizontal)
if (.not.associated(model%data_hz(link%target%read_indices%pointers(1)%p)%p)) n = n+1
if (index_column_depth==-1 .and. associated(link%target%standard_variables%first)) then
if (link%target%standard_variables%first%p%compare(standard_variables%bottom_depth)) index_column_depth = n
end if
case (domain_scalar)
if (.not.associated(model%data_scalar(link%target%read_indices%pointers(1)%p)%p)) n = n+1
end select
......@@ -48,6 +57,8 @@
link => link%next
end do
if (index_column_depth==-1) n = n + 1
! Allocate arrays to hold information on environment
allocate(environment_names(n))
allocate(environment_units(n))
......@@ -95,6 +106,13 @@
end if
link => link%next
end do
if (index_column_depth==-1) then
n = n + 1
index_column_depth = n
environment_names(n) = trim(standard_variables%bottom_depth%name)
environment_units(n) = trim(standard_variables%bottom_depth%units)
end if
end subroutine get_environment_metadata
subroutine get_couplings(model,link_list)
......
......@@ -39,6 +39,8 @@
class (type_model),private,pointer,save :: model => null()
real(8),dimension(:),pointer :: state
character(len=1024),dimension(:),allocatable :: environment_names,environment_units
integer :: index_column_depth
real(c_double),pointer :: column_depth
type (type_link_list),save :: coupling_link_list
type (type_property_dictionary),save,private :: forced_parameters,forced_couplings
......@@ -119,7 +121,8 @@
call fabm_set_domain(model)
! Retrieve arrays to hold values for environmental variables and corresponding metadata.
call get_environment_metadata(model,environment_names,environment_units)
call get_environment_metadata(model,environment_names,environment_units,index_column_depth)
column_depth => null()
call get_couplings(model,coupling_link_list)
end subroutine initialize
......@@ -171,7 +174,8 @@
call fabm_set_domain(model)
! Retrieve arrays to hold values for environmental variables and corresponding metadata.
call get_environment_metadata(model,environment_names,environment_units)
call get_environment_metadata(model,environment_names,environment_units,index_column_depth)
column_depth => null()
call get_couplings(model,coupling_link_list)
end subroutine reinitialize
......@@ -352,6 +356,7 @@
call fabm_link_bulk_data(model,environment_names(index),value)
call fabm_link_horizontal_data(model,environment_names(index),value)
call fabm_link_scalar_data(model,environment_names(index),value)
if (index==index_column_depth) column_depth => value
end subroutine link_dependency_data
subroutine link_bulk_state_data(index,value) bind(c)
......@@ -388,15 +393,20 @@
real(c_double),pointer :: pelagic_rates(:)
real(rk) :: ext
if (.not.associated(column_depth)) call driver%fatal_error('get_rates', &
'Value for environmental dependency '//trim(environment_names(index_column_depth))// &
' must be provided before calling get_rates.')
call fabm_get_light_extinction(model,ext)
call fabm_get_light(model)
call c_f_pointer(c_loc(pelagic_rates_),pelagic_rates, &
(/size(model%state_variables)+size(model%surface_state_variables)+size(model%bottom_state_variables)/))
pelagic_rates = 0.0_rk
call fabm_do_bottom(model,pelagic_rates(1:size(model%state_variables)), &
pelagic_rates(size(model%state_variables)+size(model%surface_state_variables)+1:))
call fabm_do_surface(model,pelagic_rates(1:size(model%state_variables)), &
pelagic_rates(size(model%state_variables)+1:size(model%state_variables)+size(model%surface_state_variables)))
call fabm_do_bottom(model,pelagic_rates(1:size(model%state_variables)), &
pelagic_rates(size(model%state_variables)+size(model%surface_state_variables)+1:))
pelagic_rates(1:size(model%state_variables)) = pelagic_rates(1:size(model%state_variables))/column_depth
call fabm_do(model,pelagic_rates(1:size(model%state_variables)))
! Compute rate of change in conserved quantities
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment