biodensity_feedback: initial commit

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
 

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

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

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
+!-----------------------------------------------------------------------