From c4cb78ce8e9e1c8ee1bdc614f0eb0f6a746963d7 Mon Sep 17 00:00:00 2001 From: Knut Date: Tue, 25 Aug 2020 23:53:43 +0200 Subject: [PATCH 1/4] biodensity_feedback: initial commit --- src/fabm/gotm_fabm.F90 | 24 ++++++++++++ src/gotm/gotm.F90 | 3 ++ src/meanflow/CMakeLists.txt | 1 + src/meanflow/density_correction.F90 | 57 +++++++++++++++++++++++++++++ 4 files changed, 85 insertions(+) create mode 100644 src/meanflow/density_correction.F90 diff --git a/src/fabm/gotm_fabm.F90 b/src/fabm/gotm_fabm.F90 index 1295754f3..fde0010cc 100644 --- a/src/fabm/gotm_fabm.F90 +++ b/src/fabm/gotm_fabm.F90 @@ -67,6 +67,9 @@ module gotm_fabm REALTYPE,allocatable,dimension(:,:),public :: cc_diag REALTYPE,allocatable,dimension(:), public :: cc_diag_hz +! Optional feedback to host (allocated by FABM) + REALTYPE, pointer, public :: fabm_rho_corr(:)=>NULL() ! density correction + ! Arrays for observations, relaxation times and FABM variable identifiers associated with the observations. type type_1d_state_info REALTYPE, pointer, dimension(:) :: obs => null() @@ -94,6 +97,8 @@ module gotm_fabm type (type_bulk_variable_id), save :: temp_id,salt_id,rho_id,h_id,swr_id,par_id,pres_id type (type_horizontal_variable_id),save :: lon_id,lat_id,windspeed_id,par_sf_id,cloud_id,taub_id,swr_sf_id + type (type_bulk_variable_id), save :: rho_corr_id + ! Variables to hold time spent on advection, diffusion, sink/source terms. integer(8) :: clock_adv,clock_diff,clock_source ! @@ -105,6 +110,7 @@ module gotm_fabm bioshade_feedback,bioalbedo_feedback,biodrag_feedback, & freshwater_impact,salinity_relaxation_to_freshwater_flux, & save_inputs, no_surface + logical :: biodensity_feedback=.false. ! Arrays for work, vertical movement, and cross-boundary fluxes REALTYPE,allocatable,dimension(:,:) :: ws @@ -185,6 +191,7 @@ subroutine configure_gotm_fabm_from_nml(namlst, fname) logical :: no_precipitation_dilution namelist /gotm_fabm_nml/ fabm_calc, & cnpar,w_adv_discr,ode_method,split_factor, & + biodensity_feedback, & bioshade_feedback,bioalbedo_feedback,biodrag_feedback, & repair_state,no_precipitation_dilution, & salinity_relaxation_to_freshwater_flux,save_inputs, & @@ -269,6 +276,8 @@ subroutine configure_gotm_fabm(cfg) call cfg%get(freshwater_impact, 'freshwater_impact', 'enable dilution/concentration by precipitation/evaporation', & default=.true.) ! disable to check mass conservation branch => cfg%get_child('feedbacks', 'feedbacks to physics') + call branch%get(biodensity_feedback, 'density', 'interior density correction', & + default=biodensity_feedback) call branch%get(bioshade_feedback, 'shade', 'interior light absorption', & default=.false.) call branch%get(bioalbedo_feedback, 'albedo', 'surface albedo', & @@ -485,6 +494,17 @@ subroutine init_gotm_fabm(nlev,dt,field_manager) cloud_id = model%get_horizontal_variable_id(standard_variables%cloud_area_fraction) taub_id = model%get_horizontal_variable_id(standard_variables%bottom_stress) + if (biodensity_feedback) then +! check whether used FABM version provides density_correction + rho_corr_id = model%get_bulk_variable_id("density_correction") +#if _FABM_API_VERSION_ > 0 + if (associated(rho_corr_id%variable)) then + !call model%require_data(standard_variables%density_correction) + call model%require_data(type_interior_standard_variable(name="density_correction", units="kg m-3", aggregate_variable=.true.)) + end if +#endif + end if + ! Initialize optional forcings to "off" fabm_airp => null() fabm_julianday => null() @@ -1000,6 +1020,10 @@ subroutine start_gotm_fabm(nlev, field_manager) if (total0(i) == huge(_ZERO_)) total0(i) = total(i) end do + if (biodensity_feedback .and. associated(rho_corr_id%variable)) then + fabm_rho_corr => model%get_data(rho_corr_id) + end if + end subroutine start_gotm_fabm !EOC diff --git a/src/gotm/gotm.F90 b/src/gotm/gotm.F90 index d9fce05fe..e6eb41f25 100644 --- a/src/gotm/gotm.F90 +++ b/src/gotm/gotm.F90 @@ -92,6 +92,7 @@ module gotm use gotm_fabm,only:configure_gotm_fabm,configure_gotm_fabm_from_nml,gotm_fabm_create_model,init_gotm_fabm,init_gotm_fabm_state,start_gotm_fabm,set_env_gotm_fabm,do_gotm_fabm,clean_gotm_fabm,fabm_calc use gotm_fabm,only:model_fabm=>model,standard_variables_fabm=>standard_variables use gotm_fabm, only: fabm_airp, fabm_calendar_date, fabm_julianday + use gotm_fabm, only: fabm_rho_corr use gotm_fabm_input,only: configure_gotm_fabm_input, configure_gotm_fabm_input_from_nml, init_gotm_fabm_input #endif @@ -571,6 +572,7 @@ subroutine init_gotm() #ifdef _FABM_ ! Accept the current biogeochemical state and used it to compute derived diagnostics. if (fabm_calc) call start_gotm_fabm(nlev, fm) + if (associated(fabm_rho_corr)) call density_correction(nlev, fabm_rho_corr) #endif if (list_fields) call fm%list() @@ -757,6 +759,7 @@ subroutine time_loop() #endif #ifdef _FABM_ call do_gotm_fabm(nlev,real(n,kind(_ONE_))) + if (associated(fabm_rho_corr)) call density_correction(nlev, fabm_rho_corr) #endif ! compute turbulent mixing diff --git a/src/meanflow/CMakeLists.txt b/src/meanflow/CMakeLists.txt index 482d62ad9..44238e985 100644 --- a/src/meanflow/CMakeLists.txt +++ b/src/meanflow/CMakeLists.txt @@ -2,6 +2,7 @@ add_library(meanflow buoyancy.F90 convectiveadjustment.F90 coriolis.F90 + density_correction.F90 extpressure.F90 friction.F90 intpressure.F90 diff --git a/src/meanflow/density_correction.F90 b/src/meanflow/density_correction.F90 new file mode 100644 index 000000000..08d25eb43 --- /dev/null +++ b/src/meanflow/density_correction.F90 @@ -0,0 +1,57 @@ +#include"cppdefs.h" +!----------------------------------------------------------------------- +!BOP +! +! !ROUTINE: Apply density correction +! +! !INTERFACE: + subroutine density_correction(nlev, rho_corr) +! +! !DESCRIPTION: +! +! !USES: + use meanflow, only: h, rho, buoy, NN + use meanflow, only: gravity, rho_0 + IMPLICIT NONE +! +! !INPUT PARAMETERS: + +! number of vertical layers + integer, intent(in) :: nlev + +! density correction + REALTYPE, intent(in) :: rho_corr(1:nlev) + +! +! !OUTPUT PARAMETERS: + +! +! !REVISION HISTORY: +! Original author(s): Knut Klingbeil +! +!EOP +! +! !LOCAL VARIABLES: + integer :: i + REALTYPE :: buoy_corr(1:nlev), r2b, dz +! +!----------------------------------------------------------------------- +!BOC + r2b = -gravity / rho_0 + do i=1,nlev + rho (i) = rho(i) + rho_corr(i) + buoy_corr(i) = rho_corr(i) * r2b + buoy (i) = buoy(i) + buoy_corr(i) + end do + do i=1,nlev-1 + dz = _HALF_ * ( h(i) + h(i+1) ) + NN(i) = NN(i) + ( buoy_corr(i+1) - buoy_corr(i) ) / dz + end do + + return + end subroutine density_correction +!EOC + +!----------------------------------------------------------------------- +! Copyright by the GOTM-team under the GNU Public License - www.gnu.org +!----------------------------------------------------------------------- From 935cc1f549912f8f3c878a3d66901b5790e3d196 Mon Sep 17 00:00:00 2001 From: Knut Date: Wed, 26 Aug 2020 12:58:36 +0200 Subject: [PATCH 2/4] density_feedback: turn automatic array into allocatable --- src/meanflow/density_correction.F90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/meanflow/density_correction.F90 b/src/meanflow/density_correction.F90 index 08d25eb43..1bcb370c8 100644 --- a/src/meanflow/density_correction.F90 +++ b/src/meanflow/density_correction.F90 @@ -32,11 +32,19 @@ subroutine density_correction(nlev, rho_corr) !EOP ! ! !LOCAL VARIABLES: - integer :: i - REALTYPE :: buoy_corr(1:nlev), r2b, dz + integer :: i + REALTYPE :: r2b, dz + REALTYPE, allocatable, save :: buoy_corr(:) + logical, save :: first=.true. ! !----------------------------------------------------------------------- !BOC + + if (first) then + allocate(buoy_corr(1:nlev)) + first = .false. + end if + r2b = -gravity / rho_0 do i=1,nlev rho (i) = rho(i) + rho_corr(i) From 10e76890f4752ef6365b3122d62b951ee0dc1a47 Mon Sep 17 00:00:00 2001 From: Knut Date: Wed, 26 Aug 2020 13:22:21 +0200 Subject: [PATCH 3/4] biodensity_feedback: warning for switching off biodensity_feedback --- src/fabm/gotm_fabm.F90 | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/fabm/gotm_fabm.F90 b/src/fabm/gotm_fabm.F90 index fde0010cc..bbb9eed7e 100644 --- a/src/fabm/gotm_fabm.F90 +++ b/src/fabm/gotm_fabm.F90 @@ -497,12 +497,21 @@ subroutine init_gotm_fabm(nlev,dt,field_manager) if (biodensity_feedback) then ! check whether used FABM version provides density_correction rho_corr_id = model%get_bulk_variable_id("density_correction") -#if _FABM_API_VERSION_ > 0 if (associated(rho_corr_id%variable)) then +#if _FABM_API_VERSION_ > 0 !call model%require_data(standard_variables%density_correction) call model%require_data(type_interior_standard_variable(name="density_correction", units="kg m-3", aggregate_variable=.true.)) - end if +#else +! Note: For old FABM get_data() returns null pointer. +! Maybe because require_data() was not called, but this then +! must be done BEFORE fabm_initialize(), called by default +! during create_model(). +! Note that new FABM accepts the call to require_data() +! only AFTER initialize()! + STDERR "WARNING: Density feedback from used FABM version not supported! Reset biodensity_feedback=F ..." + biodensity_feedback = .false. #endif + end if end if ! Initialize optional forcings to "off" From 8f7d4655c532f075fd3b83b1331f0e59fcb0183d Mon Sep 17 00:00:00 2001 From: Knut Date: Wed, 26 Aug 2020 15:40:20 +0200 Subject: [PATCH 4/4] biodensity_feedback: add additional check (not needed for now) --- src/fabm/gotm_fabm.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/fabm/gotm_fabm.F90 b/src/fabm/gotm_fabm.F90 index bbb9eed7e..8068faa68 100644 --- a/src/fabm/gotm_fabm.F90 +++ b/src/fabm/gotm_fabm.F90 @@ -1031,6 +1031,10 @@ subroutine start_gotm_fabm(nlev, field_manager) if (biodensity_feedback .and. associated(rho_corr_id%variable)) then fabm_rho_corr => model%get_data(rho_corr_id) + !if (associated(fabm_rho_corr,model%get_data(model%get_bulk_variable_id('zero')))) then + ! LEVEL2 "biodensity_feedback: no contribution to density_correction" + ! nullify(fabm_rho_corr) + !end if end if end subroutine start_gotm_fabm