diff --git a/src/fabm/gotm_fabm.F90 b/src/fabm/gotm_fabm.F90 index 1295754f3..8068faa68 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,26 @@ 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 (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.)) +#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" fabm_airp => null() fabm_julianday => null() @@ -1000,6 +1029,14 @@ 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) + !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 !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..1bcb370c8 --- /dev/null +++ b/src/meanflow/density_correction.F90 @@ -0,0 +1,65 @@ +#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 :: 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) + 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 +!-----------------------------------------------------------------------