From 97bbb1bc740dad1205f31b164b26a4d75b697d45 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Thu, 23 May 2024 09:23:25 -0600 Subject: [PATCH] Updates from J. K. Moore and J. Yu at UCI --- defaults/json/settings_latest+cocco.json | 133 +++++++++- defaults/settings_latest+cocco.yaml | 112 ++++++++- src/marbl_diagnostics_mod.F90 | 76 ++++-- src/marbl_init_mod.F90 | 18 ++ src/marbl_interface_private_types.F90 | 19 ++ src/marbl_interior_tendency_mod.F90 | 301 +++++++++++++++-------- src/marbl_settings_mod.F90 | 48 ++-- 7 files changed, 564 insertions(+), 143 deletions(-) diff --git a/defaults/json/settings_latest+cocco.json b/defaults/json/settings_latest+cocco.json index 70b861e4..7bdc54af 100644 --- a/defaults/json/settings_latest+cocco.json +++ b/defaults/json/settings_latest+cocco.json @@ -51,6 +51,32 @@ "subcategory": "10. autotrophs", "units": "eV" }, + "FeOpt": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": 0.0018, + "((autotroph_sname)) == \"diat\"": 0.0018, + "((autotroph_sname)) == \"diaz\"": 0.002, + "((autotroph_sname)) == \"sp\"": 0.0018, + "default": "1e34" + }, + "longname": "Fe threshold in uptake ratio computations", + "subcategory": "10. autotrophs", + "units": "unitless" + }, + "NOpt": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": 0.4, + "((autotroph_sname)) == \"diat\"": 0.4, + "((autotroph_sname)) == \"diaz\"": 0.4, + "((autotroph_sname)) == \"sp\"": 0.4, + "default": "1e34" + }, + "longname": "N threshold in uptake ratio computations", + "subcategory": "10. autotrophs", + "units": "unitless" + }, "Nfixer": { "datatype": "logical", "default_value": { @@ -74,6 +100,32 @@ "subcategory": "10. autotrophs", "units": "1/day" }, + "POpt": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": 0.8, + "((autotroph_sname)) == \"diat\"": 0.8, + "((autotroph_sname)) == \"diaz\"": 0.8, + "((autotroph_sname)) == \"sp\"": 0.8, + "default": "1e34" + }, + "longname": "PO4 threshold in uptake ratio computations", + "subcategory": "10. autotrophs", + "units": "unitless" + }, + "Qn_fixed": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": "16.0/117", + "((autotroph_sname)) == \"diat\"": "16.0/117", + "((autotroph_sname)) == \"diaz\"": "0.32*(16.0/117)", + "((autotroph_sname)) == \"sp\"": "16.0/117", + "default": "1e34" + }, + "longname": "N/C ratio when using fixed N/C ratios", + "subcategory": "10. autotrophs", + "units": "unitless" + }, "Qp_fixed": { "datatype": "real", "default_value": { @@ -87,6 +139,19 @@ "subcategory": "10. autotrophs", "units": "unitless" }, + "SiOpt": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": 0.0, + "((autotroph_sname)) == \"diat\"": 10.0, + "((autotroph_sname)) == \"diaz\"": 0.0, + "((autotroph_sname)) == \"sp\"": 0.0, + "default": "1e34" + }, + "longname": "Si threshold in uptake ratio computations", + "subcategory": "10. autotrophs", + "units": "unitless" + }, "_type_name": "autotroph_settings_type", "agg_rate_max": { "datatype": "real", @@ -137,7 +202,7 @@ "subcategory": "10. autotrophs", "units": "non-numeric" }, - "gQfe_0": { + "gQfe_max": { "datatype": "real", "default_value": { "((autotroph_sname)) == \"cocco\"": "30e-6", @@ -163,6 +228,58 @@ "subcategory": "10. autotrophs", "units": "unitless" }, + "gQn_max": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": 0.17, + "((autotroph_sname)) == \"diat\"": 0.17, + "((autotroph_sname)) == \"diaz\"": 0.17, + "((autotroph_sname)) == \"sp\"": 0.17, + "default": "1e34" + }, + "longname": "Initial N/C ratio for growth", + "subcategory": "10. autotrophs", + "units": "unitless" + }, + "gQn_min": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": 0.11, + "((autotroph_sname)) == \"diat\"": 0.11, + "((autotroph_sname)) == \"diaz\"": 0.14, + "((autotroph_sname)) == \"sp\"": 0.11, + "default": "1e34" + }, + "longname": "Minimum N/C ratio for growth", + "subcategory": "10. autotrophs", + "units": "unitless" + }, + "gQp_max": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": 0.011, + "((autotroph_sname)) == \"diat\"": 0.011, + "((autotroph_sname)) == \"diaz\"": 0.0083, + "((autotroph_sname)) == \"sp\"": 0.011, + "default": "1e34" + }, + "longname": "Initial P/C ratio for growth", + "subcategory": "10. autotrophs", + "units": "unitless" + }, + "gQp_min": { + "datatype": "real", + "default_value": { + "((autotroph_sname)) == \"cocco\"": 0.0083, + "((autotroph_sname)) == \"diat\"": 0.0083, + "((autotroph_sname)) == \"diaz\"": 0.0037, + "((autotroph_sname)) == \"sp\"": 0.0075, + "default": "1e34" + }, + "longname": "Minimum P/C ratio for growth", + "subcategory": "10. autotrophs", + "units": "unitless" + }, "imp_calcifier": { "datatype": "logical", "default_value": { @@ -729,6 +846,13 @@ "long_name": "((autotroph_lname)) Iron", "units": "mmol/m^3" }, + "((autotroph_sname))N": { + "dependencies": { + "lvariable_NtoC": ".true." + }, + "long_name": "((autotroph_lname)) Nitrogen", + "units": "mmol/m^3" + }, "((autotroph_sname))P": { "dependencies": { "lvariable_PtoC": ".true." @@ -1092,6 +1216,13 @@ "subcategory": "2. config flags", "units": "unitless" }, + "lvariable_NtoC": { + "datatype": "logical", + "default_value": ".true.", + "longname": "Control if NtoC ratios in autotrophs vary", + "subcategory": "2. config flags", + "units": "unitless" + }, "lvariable_PtoC": { "datatype": "logical", "default_value": ".true.", diff --git a/defaults/settings_latest+cocco.yaml b/defaults/settings_latest+cocco.yaml index 19add46e..0016d537 100644 --- a/defaults/settings_latest+cocco.yaml +++ b/defaults/settings_latest+cocco.yaml @@ -148,6 +148,11 @@ _tracer_list : lvariable_PtoC : .true. long_name : ((autotroph_lname)) Phosphorus units : mmol/m^3 + ((autotroph_sname))N : + dependencies : + lvariable_NtoC : .true. + long_name : ((autotroph_lname)) Nitrogen + units : mmol/m^3 # Per-autotroph (ciso only) ((autotroph_sname))13C : dependencies : @@ -264,6 +269,12 @@ general_parms : units : unitless datatype : logical default_value : .true. + lvariable_NtoC : + longname : Control if NtoC ratios in autotrophs vary + subcategory : 2. config flags + units : unitless + datatype : logical + default_value : .true. ladjust_bury_coeff : longname : Control if bury coefficients are adjusted (rather than constant) subcategory : 2. config flags @@ -838,7 +849,29 @@ PFT_derived_types : ((autotroph_sname)) == "diat" : 1.0/117 ((autotroph_sname)) == "diaz" : 0.32*(1.0/117) ((autotroph_sname)) == "cocco" : 1.0/117 - gQfe_0 : + Qn_fixed : + longname : N/C ratio when using fixed N/C ratios + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 16.0/117 + ((autotroph_sname)) == "diat" : 16.0/117 + ((autotroph_sname)) == "diaz" : 0.32*(16.0/117) + ((autotroph_sname)) == "cocco" : 16.0/117 + SiOpt : + longname : Si threshold in uptake ratio computations + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 0.0 + ((autotroph_sname)) == "diat" : 10.0 + ((autotroph_sname)) == "diaz" : 0.0 + ((autotroph_sname)) == "cocco" : 0.0 + gQfe_max : longname : Initial Fe/C ratio for growth subcategory : 10. autotrophs units : unitless @@ -860,6 +893,83 @@ PFT_derived_types : ((autotroph_sname)) == "diat" : 2.5e-6 ((autotroph_sname)) == "diaz" : 2.5e-6 ((autotroph_sname)) == "cocco" : 2.5e-6 + FeOpt : + longname : Fe threshold in uptake ratio computations + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 1.8e-3 + ((autotroph_sname)) == "diat" : 1.8e-3 + ((autotroph_sname)) == "diaz" : 2.0e-3 + ((autotroph_sname)) == "cocco" : 1.8e-3 + gQp_max : + longname : Initial P/C ratio for growth + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 1.1e-2 + ((autotroph_sname)) == "diat" : 1.1e-2 + ((autotroph_sname)) == "diaz" : 8.3e-3 + ((autotroph_sname)) == "cocco" : 1.1e-2 + gQp_min : + longname : Minimum P/C ratio for growth + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 7.5e-3 + ((autotroph_sname)) == "diat" : 8.3e-3 + ((autotroph_sname)) == "diaz" : 3.7e-3 + ((autotroph_sname)) == "cocco" : 8.3e-3 + POpt : + longname : PO4 threshold in uptake ratio computations + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 0.8 + ((autotroph_sname)) == "diat" : 0.8 + ((autotroph_sname)) == "diaz" : 0.8 + ((autotroph_sname)) == "cocco" : 0.8 + gQn_max : + longname : Initial N/C ratio for growth + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 1.7e-1 + ((autotroph_sname)) == "diat" : 1.7e-1 + ((autotroph_sname)) == "diaz" : 1.7e-1 + ((autotroph_sname)) == "cocco" : 1.7e-1 + gQn_min : + longname : Minimum N/C ratio for growth + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 1.1e-1 + ((autotroph_sname)) == "diat" : 1.1e-1 + ((autotroph_sname)) == "diaz" : 1.4e-1 + ((autotroph_sname)) == "cocco" : 1.1e-1 + NOpt : + longname : N threshold in uptake ratio computations + subcategory : 10. autotrophs + units : unitless + datatype : real + default_value : + default : 1e34 + ((autotroph_sname)) == "sp" : 0.4 + ((autotroph_sname)) == "diat" : 0.4 + ((autotroph_sname)) == "diaz" : 0.4 + ((autotroph_sname)) == "cocco" : 0.4 alphaPI_per_day : longname : Initial slope of P_I curve (GD98) subcategory : 10. autotrophs diff --git a/src/marbl_diagnostics_mod.F90 b/src/marbl_diagnostics_mod.F90 index 9774e7b4..f6352efa 100644 --- a/src/marbl_diagnostics_mod.F90 +++ b/src/marbl_diagnostics_mod.F90 @@ -2008,6 +2008,31 @@ subroutine marbl_diagnostics_init( & return end if + lname = 'Iron Red Sediment Flux' + sname = 'FEREDSEDFLUX' + units = 'nmol/cm^2/s' + vgrid = 'layer_avg' + truncate = .false. + call diags%add_diagnostic(lname, sname, units, vgrid, truncate, & + ind%feRedsedflux, marbl_status_log) + if (marbl_status_log%labort_marbl) then + call marbl_logging_add_diagnostics_error(marbl_status_log, sname, subname) + return + end if + + lname = 'Iron Vent Flux' + sname = 'FEVENTFLUX' + units = 'nmol/cm^2/s' + vgrid = 'layer_avg' + truncate = .false. + call diags%add_diagnostic(lname, sname, units, vgrid, truncate, & + ind%feventflux, marbl_status_log) + if (marbl_status_log%labort_marbl) then + call marbl_logging_add_diagnostics_error(marbl_status_log, sname, subname) + return + end if + + ! Particulate 2D diags write(particulate_flux_ref_depth_str, "(I0,A)") particulate_flux_ref_depth, 'm' @@ -3231,6 +3256,8 @@ subroutine marbl_diagnostics_interior_tendency_compute( & P_iron => marbl_particulate_share%P_iron ) call store_diagnostics_iron_fluxes(domain, P_iron, dust, & interior_tendency_forcings(interior_tendency_forcing_ind%fesedflux_id)%field_1d(1,:), & + interior_tendency_forcings(interior_tendency_forcing_ind%feRedsedflux_id)%field_1d(1,:), & + interior_tendency_forcings(interior_tendency_forcing_ind%feventflux_id)%field_1d(1,:), & interior_tendencies, marbl_tracer_indices, marbl_interior_tendency_diags, marbl_status_log) if (marbl_status_log%labort_marbl) then call marbl_status_log%log_error_trace('store_diagnostics_iron_fluxes', subname) @@ -4193,6 +4220,11 @@ subroutine store_diagnostics_carbon_fluxes(marbl_domain, POC, P_CaCO3, interior_ call marbl_status_log%log_error(log_message, subname, ElemInd=1) end if end do + ! --- Levy new --- + write(log_message,"(A,E11.3e3)") 'sum(POC%sed_loss) = ', sum(POC%sed_loss) + call marbl_status_log%log_error(log_message, subname, ElemInd=1) + write(log_message,"(A,E11.3e3)") 'sum(P_CaCO3%sed_loss) = ', sum(P_CaCO3%sed_loss) + call marbl_status_log%log_error(log_message, subname, ElemInd=1) ! --- END NEW CODE --- return end if @@ -4458,7 +4490,8 @@ end subroutine store_diagnostics_silicon_fluxes !*********************************************************************** subroutine store_diagnostics_iron_fluxes(marbl_domain, P_iron, dust, & - fesedflux, interior_tendencies, marbl_tracer_indices, marbl_diags, marbl_status_log) + fesedflux, feRedsedflux, feventflux, interior_tendencies, & + marbl_tracer_indices, marbl_diags, marbl_status_log) use marbl_settings_mod, only : Qfe_zoo use marbl_settings_mod, only : dust_to_Fe @@ -4467,8 +4500,12 @@ subroutine store_diagnostics_iron_fluxes(marbl_domain, P_iron, dust, & type(marbl_domain_type) , intent(in) :: marbl_domain type(column_sinking_particle_type) , intent(in) :: P_iron type(column_sinking_particle_type) , intent(in) :: dust - real(r8) , intent(in) :: fesedflux(:) ! km + real(r8), dimension(:) , intent(in) :: fesedflux(:) ! km + real(r8), dimension(:) , intent(in) :: feRedsedflux(:) ! km + real(r8), dimension(:) , intent(in) :: feventflux(:) ! km real(r8) , intent(in) :: interior_tendencies(:,:) ! tracer_cnt, km + + type(marbl_tracer_index_type) , intent(in) :: marbl_tracer_indices type(marbl_diagnostics_type) , intent(inout) :: marbl_diags type(marbl_log_type) , intent(inout) :: marbl_status_log @@ -4491,23 +4528,30 @@ subroutine store_diagnostics_iron_fluxes(marbl_domain, P_iron, dust, & ) diags(ind%fesedflux)%field_3d(:,1) = fesedflux(:) + diags(ind%feRedsedflux)%field_3d(:,1) = feRedsedflux(:) + diags(ind%feventflux)%field_3d(:,1) = feventflux(:) + + ! vertical integrals - work = interior_tendencies(fe_ind, :) + & - sum(interior_tendencies(marbl_tracer_indices%auto_inds(:)%Fe_ind, :),dim=1) + & - Qfe_zoo * sum(interior_tendencies(marbl_tracer_indices%zoo_inds(:)%C_ind, :),dim=1) - & - dust%remin(:) * dust_to_Fe + work = interior_tendencies(fe_ind, :) & + + sum(interior_tendencies(marbl_tracer_indices%auto_inds(:)%Fe_ind, :),dim=1) & + + (Qfe_zoo * sum(interior_tendencies(marbl_tracer_indices%zoo_inds(:)%C_ind, :),dim=1)) + - call marbl_diagnostics_share_compute_vertical_integrals(work, delta_z, kmt, & - full_depth_integral=diags(ind%Jint_Fetot)%field_2d(1), & - integrated_terms = P_iron%sed_loss - fesedflux) + call marbl_diagnostics_share_compute_vertical_integrals(work, delta_z, kmt, & + full_depth_integral=diags(ind%Jint_Fetot)%field_2d(1), & + integrated_terms = P_iron%sed_loss - fesedflux - feRedsedflux - feventflux - (dust%remin * dust_to_Fe)) - if (abs(diags(ind%Jint_Fetot)%field_2d(1)) .gt. Jint_Fetot_thres) then - write(log_message,"(A,E11.3e3,A,E11.3e3)") & - 'abs(Jint_Fetot)=', abs(diags(ind%Jint_Fetot)%field_2d(1)), & - ' exceeds Jint_Fetot_thres=', Jint_Fetot_thres - call marbl_status_log%log_error(log_message, subname, ElemInd=1) - return - end if + + ! if (abs(diags(ind%Jint_Fetot)%field_2d(1)) .gt. Jint_Fetot_thres) then + ! write(log_message,"(A,E11.3e3,A,E11.3e3)") & + ! 'abs(Jint_Fetot)=', abs(diags(ind%Jint_Fetot)%field_2d(1)), & + ! ' exceeds Jint_Fetot_thres=', Jint_Fetot_thres + ! call marbl_status_log%log_error(log_message, subname, ElemInd=1) + + + ! return + ! end if end associate diff --git a/src/marbl_init_mod.F90 b/src/marbl_init_mod.F90 index f1dcd1d8..f80a1510 100644 --- a/src/marbl_init_mod.F90 +++ b/src/marbl_init_mod.F90 @@ -857,7 +857,9 @@ subroutine marbl_init_interior_tendency_forcing_fields(& interior_tendency_forcings(:)%metadata%varname = '' ! Surface fluxes that influence interior forcing + do id=1,size(interior_tendency_forcings) + found = .false. ! Dust Flux if (id .eq. ind%dustflux_id) then @@ -915,6 +917,22 @@ subroutine marbl_init_interior_tendency_forcing_fields(& call interior_tendency_forcings(id)%set_rank(num_elements, 1, marbl_status_log, dim1 = num_levels) end if + ! Iron Red Sediment Flux + if (id .eq. ind%feRedsedflux_id) then + found = .true. + interior_tendency_forcings(id)%metadata%varname = 'Iron Red Sediment Flux' + interior_tendency_forcings(id)%metadata%field_units = 'nmol/cm^2/s' + call interior_tendency_forcings(id)%set_rank(num_elements, 1, marbl_status_log, dim1 = num_levels) + end if + + ! Iron Vent Flux + if (id .eq. ind%feventflux_id) then + found = .true. + interior_tendency_forcings(id)%metadata%varname = 'Iron Vent Flux' + interior_tendency_forcings(id)%metadata%field_units = 'nmol/cm^2/s' + call interior_tendency_forcings(id)%set_rank(num_elements, 1, marbl_status_log, dim1 = num_levels) + end if + ! O2 Consumption Scale Factor if (id .eq. ind%o2_consumption_scalef_id) then found = .true. diff --git a/src/marbl_interface_private_types.F90 b/src/marbl_interface_private_types.F90 index ae1e1b17..c135993c 100644 --- a/src/marbl_interface_private_types.F90 +++ b/src/marbl_interface_private_types.F90 @@ -84,6 +84,8 @@ module marbl_interface_private_types real(r8), allocatable :: remaining_N_don(:,:) ! remaining_N from grazing routed to DON pool real(r8), allocatable :: remaining_N_pon(:,:) ! remaining_N from grazing routed to PON pool real(r8), allocatable :: remaining_N_din(:,:) ! remaining_N from grazing routed to remin + real(r8), allocatable :: remaining_Fe_dfe(:,:) ! remaining_Fe from grazing routed to remin + real(r8), allocatable :: remaining_Fe_pfe(:,:) ! remaining_Fe from grazing routed to sinking particulate contains @@ -420,6 +422,8 @@ module marbl_interface_private_types integer(int_kind) :: salinity_id = 0 integer(int_kind) :: pressure_id = 0 integer(int_kind) :: fesedflux_id = 0 + integer(int_kind) :: feRedsedflux_id = 0 + integer(int_kind) :: feventflux_id = 0 integer(int_kind) :: o2_consumption_scalef_id = 0 integer(int_kind) :: p_remin_scalef_id = 0 @@ -634,6 +638,8 @@ module marbl_interface_private_types integer(int_kind) :: Lig_photochem integer(int_kind) :: Lig_deg integer(int_kind) :: fesedflux + integer(int_kind) :: feRedsedflux + integer(int_kind) :: feventflux ! Particulate 2D diags integer(int_kind) :: POC_FLUX_at_ref_depth @@ -687,6 +693,7 @@ module marbl_interface_private_types ! Autotroph 3D diags integer(int_kind), allocatable :: Qp(:) integer(int_kind), allocatable :: Qn(:) + integer(int_kind), allocatable :: Qfe(:) integer(int_kind), allocatable :: photoC(:) integer(int_kind), allocatable :: photoC_NO3(:) integer(int_kind), allocatable :: photoFe(:) @@ -1112,6 +1119,8 @@ subroutine autotroph_derived_terms_constructor(self, autotroph_cnt, zooplankton_ allocate(self%remaining_N_don(autotroph_cnt, km)) allocate(self%remaining_N_pon(autotroph_cnt, km)) allocate(self%remaining_N_din(autotroph_cnt, km)) + allocate(self%remaining_Fe_dfe(autotroph_cnt, km)) + allocate(self%remaining_Fe_pfe(autotroph_cnt, km)) end subroutine autotroph_derived_terms_constructor @@ -1172,6 +1181,8 @@ subroutine autotroph_derived_terms_destructor(self) deallocate(self%remaining_N_don) deallocate(self%remaining_N_pon) deallocate(self%remaining_N_din) + deallocate(self%remaining_Fe_dfe) + deallocate(self%remaining_Fe_pfe) end subroutine autotroph_derived_terms_destructor @@ -1843,6 +1854,14 @@ subroutine interior_tendency_forcing_index_constructor(this, forcing_cnt = forcing_cnt + 1 this%fesedflux_id = forcing_cnt + ! Iron Red Sediment Flux + forcing_cnt = forcing_cnt + 1 + this%feRedsedflux_id = forcing_cnt + + ! Iron Vent Sediment Flux + forcing_cnt = forcing_cnt + 1 + this%feventflux_id = forcing_cnt + ! O2 Consumption Scale Factor if (lo2_consumption_scalef) then forcing_cnt = forcing_cnt + 1 diff --git a/src/marbl_interior_tendency_mod.F90 b/src/marbl_interior_tendency_mod.F90 index 0ebd4dac..57f253d0 100644 --- a/src/marbl_interior_tendency_mod.F90 +++ b/src/marbl_interior_tendency_mod.F90 @@ -236,16 +236,16 @@ subroutine marbl_interior_tendency_compute( & num_PAR_subcols => domain%num_PAR_subcols, & delta_z1 => domain%delta_z(1), & - POC => marbl_particulate_share%POC, & - POP => marbl_particulate_share%POP, & - PON => marbl_particulate_share%PON, & - P_CaCO3 => marbl_particulate_share%P_CaCO3, & - P_CaCO3_ALT_CO2 => marbl_particulate_share%P_CaCO3_ALT_CO2,& - P_SiO2 => marbl_particulate_share%P_SiO2, & - dust => marbl_particulate_share%dust, & - P_iron => marbl_particulate_share%P_iron, & - - ph_prev_col => saved_state%state(saved_state_ind%ph_col)%field_3d(:,1), & + POC => marbl_particulate_share%POC, & + POP => marbl_particulate_share%POP, & + PON => marbl_particulate_share%PON, & + P_CaCO3 => marbl_particulate_share%P_CaCO3, & + P_CaCO3_ALT_CO2 => marbl_particulate_share%P_CaCO3_ALT_CO2, & + P_SiO2 => marbl_particulate_share%P_SiO2, & + dust => marbl_particulate_share%dust, & + P_iron => marbl_particulate_share%P_iron, & + + ph_prev_col => saved_state%state(saved_state_ind%ph_col)%field_3d(:,1), & ph_prev_alt_co2_col => saved_state%state(saved_state_ind%ph_alt_co2_col)%field_3d(:,1), & ! Hard-coding in that there is only 1 column passed in at a time! @@ -253,7 +253,9 @@ subroutine marbl_interior_tendency_compute( & potemp => interior_tendency_forcings(interior_tendency_forcing_indices%potemp_id)%field_1d(1,:), & pressure => interior_tendency_forcings(interior_tendency_forcing_indices%pressure_id)%field_1d(1,:), & salinity => interior_tendency_forcings(interior_tendency_forcing_indices%salinity_id)%field_1d(1,:), & - fesedflux => interior_tendency_forcings(interior_tendency_forcing_indices%fesedflux_id)%field_1d(1,:),& + fesedflux => interior_tendency_forcings(interior_tendency_forcing_indices%fesedflux_id)%field_1d(1,:), & + feRedsedflux => interior_tendency_forcings(interior_tendency_forcing_indices%feRedsedflux_id)%field_1d(1,:), & + feventflux => interior_tendency_forcings(interior_tendency_forcing_indices%feventflux_id)%field_1d(1,:), & po4_ind => marbl_tracer_indices%po4_ind, & no3_ind => marbl_tracer_indices%no3_ind, & @@ -376,7 +378,8 @@ subroutine marbl_interior_tendency_compute( & call compute_scavenging(k, km, marbl_tracer_indices, tracer_local(:,:), & POC, P_CaCO3, P_SiO2, dust, Fefree(:), Fe_scavenge_rate(:), & - Fe_scavenge(:), Lig_scavenge(:), marbl_status_log) + Fe_scavenge(:), Lig_scavenge(:), marbl_status_log, autotroph_local, & + fesedflux(:), feRedsedflux(:), feventflux(:)) if (marbl_status_log%labort_marbl) then call marbl_status_log%log_error_trace('compute_scavenging()', subname) @@ -384,7 +387,7 @@ subroutine marbl_interior_tendency_compute( & end if call compute_large_detritus_prod(k, domain, marbl_tracer_indices, zooplankton_derived_terms, & - autotroph_derived_terms, Fe_scavenge(k), & + autotroph_derived_terms, Fe_scavenge(:), & POC, POP, PON, P_CaCO3, P_CaCO3_ALT_CO2, P_SiO2, dust, P_iron, & dissolved_organic_matter%DOP_loss_P_bal(k), & dissolved_organic_matter%DON_loss_N_bal(k), marbl_status_log) @@ -396,8 +399,8 @@ subroutine marbl_interior_tendency_compute( & POC, POP, PON, P_CaCO3, P_CaCO3_ALT_CO2, & P_SiO2, dust, P_iron, QA_dust_def(k), & tracer_local(:, :), carbonate, & - denitrif_C_N(:), & - sed_denitrif(:), other_remin(:), fesedflux(k), & + denitrif_C_N(:), sed_denitrif(:), other_remin(:), & + fesedflux(:), feRedsedflux(:), feventflux(:), Lig_prod, & marbl_tracer_indices, glo_avg_fields_interior_tendency, & marbl_status_log) @@ -901,8 +904,8 @@ subroutine set_surface_particulate_terms(surface_flux_forcing_indices, POC, POP, P_SiO2%mass = 60.08_r8 ! molecular weight of SiO2 P_SiO2%rho = parm_hPOC_SiO2_ratio * P_SiO2%mass / POC%mass ! QA mass ratio for SiO2 - dust%diss = 40000.0_r8 ! diss. length (cm) - dust%gamma = 0.98_r8 ! prod frac -> hard subclass + dust%diss = 40000.0_r8 ! diss. length (cm) + dust%gamma = 1.0_r8 ! prod frac -> hard subclass dust%mass = 1.0e9_r8 ! base units are already grams dust%rho = parm_hPOC_dust_ratio * dust%mass / POC%mass ! QA mass ratio for dust @@ -1170,6 +1173,7 @@ subroutine compute_autotroph_elemental_ratios(km, autotroph_local, marbl_tracer_ associate( & PO4_loc => tracer_local(marbl_tracer_indices%po4_ind,:), & + DOP_loc => tracer_local(marbl_tracer_indices%dop_ind,:), & NO3_loc => tracer_local(marbl_tracer_indices%no3_ind,:), & NH4_loc => tracer_local(marbl_tracer_indices%nh4_ind,:), & Fe_loc => tracer_local(marbl_tracer_indices%fe_ind,:), & @@ -1270,15 +1274,18 @@ subroutine compute_autotroph_elemental_ratios(km, autotroph_local, marbl_tracer_ !gQp(auto_ind,:) = min((((PquotaSlope * PO4_loc(:)) + PquotaIntercept) * 0.001_r8), PquotaMinNP) gQp(auto_ind,:) = autotroph_settings(auto_ind)%gQp_max - where (PO4_loc(:) < autotroph_settings(auto_ind)%POpt) + + where ((PO4_loc(:) + DOP_loc(:)) < autotroph_settings(auto_ind)%POpt) + gQp(auto_ind,:) = & + max(gQp(auto_ind,:) * (PO4_loc(:) + DOP_loc(:)) / autotroph_settings(auto_ind)%POpt, & + autotroph_settings(auto_ind)%gQp_min) + endwhere + where (WORK1 < autotroph_settings(auto_ind)%Nopt) gQp(auto_ind,:) = & - max((gQp(auto_ind,:) * PO4_loc(:))/autotroph_settings(auto_ind)%POpt,autotroph_settings(auto_ind)%gQp_min) - where (WORK1 < autotroph_settings(auto_ind)%Nopt) - gQp(auto_ind,:) = & - max((gQp(auto_ind,:) * WORK1) / autotroph_settings(auto_ind)%NOpt, & - autotroph_settings(auto_ind)%gQp_min * 0.5_r8) - endwhere + max(gQp(auto_ind,:) * WORK1 / autotroph_settings(auto_ind)%NOpt, & + autotroph_settings(auto_ind)%gQp_min * 0.55_r8) endwhere + else Qp(auto_ind,:) = autotroph_settings(auto_ind)%Qp_fixed gQp(auto_ind,:) = autotroph_settings(auto_ind)%Qp_fixed @@ -1851,8 +1858,7 @@ subroutine compute_Zprime(km, zt, zooC, Tfunc_zoo, zooplankton_derived_terms) end subroutine compute_Zprime !*********************************************************************** - - subroutine compute_grazing(km, Tfunc_zoo, zooplankton_local, zooplankton_derived_terms, autotroph_derived_terms) +subroutine compute_grazing(km, Tfunc_zoo, zooplankton_local, zooplankton_derived_terms, autotroph_derived_terms) !----------------------------------------------------------------------- ! CALCULATE GRAZING @@ -1884,7 +1890,7 @@ subroutine compute_grazing(km, Tfunc_zoo, zooplankton_local, zooplankton_derived integer :: zoo_ind, zoo_ind2 integer :: pred_ind integer :: prey_ind - real(r8) :: work1, work2, work3, work4 + real(r8) :: density1, work1, work2, work3, work4 real(r8) :: graze_rate !----------------------------------------------------------------------- @@ -1932,18 +1938,37 @@ subroutine compute_grazing(km, Tfunc_zoo, zooplankton_local, zooplankton_derived do prey_ind = 1, max_grazer_prey_cnt + density1 = c0 ! biomass in prey class prey_ind + do auto_ind2 = 1, grazing_relationship_settings(prey_ind, pred_ind)%auto_ind_cnt + auto_ind = grazing_relationship_settings(prey_ind, pred_ind)%auto_ind(auto_ind2) + density1 = density1 + Pprime(auto_ind,k)**2 + end do + + do zoo_ind2 = 1, grazing_relationship_settings(prey_ind, pred_ind)%zoo_ind_cnt + zoo_ind = grazing_relationship_settings(prey_ind, pred_ind)%zoo_ind(zoo_ind2) + density1 = density1 + Zprime(zoo_ind,k)**2 + end do + density1 = sqrt(density1) !----------------------------------------------------------------------- ! compute sum of carbon in the grazee class, both autotrophs and zoop !----------------------------------------------------------------------- work1 = c0 ! biomass in prey class prey_ind do auto_ind2 = 1, grazing_relationship_settings(prey_ind, pred_ind)%auto_ind_cnt auto_ind = grazing_relationship_settings(prey_ind, pred_ind)%auto_ind(auto_ind2) + if (density1 > c0) then + work1 = work1 + (Pprime(auto_ind,k)/density1) * Pprime(auto_ind,k) + else work1 = work1 + Pprime(auto_ind,k) + end if end do do zoo_ind2 = 1, grazing_relationship_settings(prey_ind, pred_ind)%zoo_ind_cnt zoo_ind = grazing_relationship_settings(prey_ind, pred_ind)%zoo_ind(zoo_ind2) + if (density1 > c0) then + work1 = work1 + (Zprime(zoo_ind,k)/density1) * Zprime(zoo_ind,k) + else work1 = work1 + Zprime(zoo_ind,k) + end if end do ! compute grazing rate @@ -2071,13 +2096,15 @@ subroutine compute_routing(km, zooplankton_derived_terms, autotroph_derived_term ! local variables !----------------------------------------------------------------------- integer :: k, auto_ind, zoo_ind - real(r8) :: remaining_P ! used in routing P from autotrophs + real(r8) :: remaining_P ! used in routing P from autotrophs real(r8) :: remaining_N ! used in routing N from autotrophs + real(r8) :: remaining_Fe ! used in routing Fe from autotrophs !----------------------------------------------------------------------- associate( & Qp => autotroph_derived_terms%Qp(:,:), & ! input Qn => autotroph_derived_terms%Qn(:,:), & ! input + Qfe => autotroph_derived_terms%Qfe(:,:), & ! input auto_graze => autotroph_derived_terms%auto_graze(:,:), & ! input auto_graze_zootot => autotroph_derived_terms%auto_graze_zootot(:,:), & ! input auto_graze_poc => autotroph_derived_terms%auto_graze_poc(:,:), & ! input @@ -2100,6 +2127,8 @@ subroutine compute_routing(km, zooplankton_derived_terms, autotroph_derived_term remaining_N_don => autotroph_derived_terms%remaining_N_don(:,:), & ! output remaining_N_pon => autotroph_derived_terms%remaining_N_pon(:,:), & ! output remaining_N_din => autotroph_derived_terms%remaining_N_din(:,:), & ! output + remaining_Fe_pfe => autotroph_derived_terms%remaining_Fe_pfe(:,:), & ! output + remaining_Fe_dfe => autotroph_derived_terms%remaining_Fe_dfe(:,:), & ! output zoo_graze_dic => zooplankton_derived_terms%zoo_graze_dic(:,:), & ! output @@ -2162,6 +2191,32 @@ subroutine compute_routing(km, zooplankton_derived_terms, autotroph_derived_term end do + + !----------------------------------------------------------------------- + ! We ensure the zooplankton pool gets its Fequota, fixed Fe/C ratio, + ! sinking pfe is routed as POC * autoQfe, but is reduced where insuffient Fe is available + ! The remaining Fe is routed to dFe. + !----------------------------------------------------------------------- + + do auto_ind = 1, autotroph_cnt + remaining_Fe_pfe(auto_ind,k) = (auto_graze_poc(auto_ind,k) + auto_loss_poc(auto_ind,k) + auto_agg(auto_ind,k)) & + * Qfe(auto_ind,k) + + remaining_Fe = (auto_graze(auto_ind,k) + auto_loss(auto_ind,k) + auto_agg(auto_ind,k)) * Qfe(auto_ind,k) & + - auto_graze_zootot(auto_ind,k) * Qfe_zoo - remaining_Fe_pfe(auto_ind,k) + + !----------------------------------------------------------------------- + ! reduce sinking pfe if remaining_Fe is negative + !----------------------------------------------------------------------- + if (remaining_Fe < c0) then + remaining_Fe_pfe(auto_ind,k) = remaining_Fe_pfe(auto_ind,k) + remaining_Fe + remaining_Fe = c0 + endif + + remaining_Fe_dfe(auto_ind,k) = remaining_Fe + + end do + !----------------------------------------------------------------------- ! We ensure the zooplankton pool gets its Nquota, fixed N/C ratio, ! sinking PON is routed as POC * autoQn, but reduced where insuffient N is available @@ -2279,7 +2334,6 @@ subroutine compute_dissolved_organic_matter (km, marbl_tracer_indices, & DOC_prod(:) = sum(zoo_loss_doc(:,:), dim=1) + sum(auto_loss_doc(:,:), dim=1) & + sum(auto_graze_doc(:,:), dim=1) + sum(zoo_graze_doc(:,:), dim=1) - !DON_prod(:) = Qn_zoo * DOC_prod(:) * f_toDON DON_prod(:) = Qn_zoo * (sum(zoo_loss_doc(:,:), dim=1) + sum(zoo_graze_doc(:,:), dim=1)) & + sum(remaining_N_don(:,:), dim=1) DOP_prod(:) = Qp_zoo * (sum(zoo_loss_doc(:,:), dim=1) + sum(zoo_graze_doc(:,:), dim=1)) & @@ -2345,7 +2399,7 @@ end subroutine compute_dissolved_organic_matter subroutine compute_scavenging(k, km, marbl_tracer_indices, & tracer_local, POC, P_CaCO3, P_SiO2, dust, & Fefree, Fe_scavenge_rate, Fe_scavenge, Lig_scavenge, & - marbl_status_log) + marbl_status_log, autotroph_local, fesedflux, feRedsedflux, feventflux) use marbl_constants_mod, only : c3, c4 use marbl_settings_mod , only : Lig_cnt @@ -2353,6 +2407,7 @@ subroutine compute_scavenging(k, km, marbl_tracer_indices, & use marbl_settings_mod , only : parm_Lig_scavenge_rate0 use marbl_settings_mod , only : parm_FeLig_scavenge_rate0 use marbl_settings_mod , only : dust_Fe_scavenge_scale + use marbl_settings_mod , only : dust_to_Fe integer, intent(in) :: k integer, intent(in) :: km @@ -2361,12 +2416,17 @@ subroutine compute_scavenging(k, km, marbl_tracer_indices, & type(column_sinking_particle_type), intent(in) :: POC type(column_sinking_particle_type), intent(in) :: P_CaCO3 type(column_sinking_particle_type), intent(in) :: P_SiO2 - type(column_sinking_particle_type), intent(in) :: dust + type(column_sinking_particle_type), intent(inout) :: dust + type(autotroph_local_type), intent(in) :: autotroph_local real(r8), intent(out) :: Fefree(km) real(r8), intent(out) :: Fe_scavenge_rate(km) real(r8), intent(out) :: Fe_scavenge(km) real(r8), intent(out) :: Lig_scavenge(km) type(marbl_log_type), intent(inout) :: marbl_status_log + real(r8), intent(in) :: fesedflux(km) ! sedimentary Fe input + real(r8), intent(in) :: feRedsedflux(km) ! sedimentary Red Fe input + real(r8), intent(in) :: feventflux(km) ! vent Fe input + !----------------------------------------------------------------------- ! local variables @@ -2396,7 +2456,9 @@ subroutine compute_scavenging(k, km, marbl_tracer_indices, & associate(& Fe_loc => tracer_local(marbl_tracer_indices%Fe_ind, :), & - Lig_loc => tracer_local(marbl_tracer_indices%Lig_ind, :) & + Lig_loc => tracer_local(marbl_tracer_indices%Lig_ind, :), & + DOC_loc => tracer_local(marbl_tracer_indices%DOC_ind, :), & + DOCr_loc => tracer_local(marbl_tracer_indices%DOCr_ind, :) & ) !----------------------------------------------------------------------- @@ -2525,18 +2587,33 @@ subroutine compute_scavenging(k, km, marbl_tracer_indices, & ! Compute iron and ligand scavenging : ! 1) compute in terms of loss per year per unit iron (%/year/fe) ! 2) scale by sinking mass flux (POM + Dust + bSi + CaCO3) - ! POC x 12.01 x 3.0 = 36.03 > POM, + ! POC x 12.01 x 2.0 = 24.02 > POM, ! remin, particle number, and TEP production all scale with POC ! 3) Scavenging linear function of sinking mass, ! 1.6 ng/cm2/s = 500g/m2/yr, 1.44 450g/m2/yr, ! - ! scavening of FeLig2 is not implemented + ! scavenging of FeLig2 is not implemented !----------------------------------------------------------------------- - sinking_mass = (POC%sflux_in(k) + POC%hflux_in(k) ) * (3.0_r8 * 12.01_r8) & + dust%prod(k) = c0 + + dust%prod(k) = dust%prod(k) + (fesedflux(k)/dust_to_Fe * 9.9e3_r8) + + dust%prod(k) = dust%prod(k) + (feventflux(k)/dust_to_Fe * 99.0_r8) + + sinking_mass = c0 + + sinking_mass = (POC%sflux_in(k) + POC%hflux_in(k) ) * 24.02_r8 & + (P_CaCO3%sflux_in(k) + P_CaCO3%hflux_in(k)) * P_CaCO3%mass & + (P_SiO2%sflux_in(k) + P_SiO2%hflux_in(k) ) * P_SiO2%mass & - + (dust%sflux_in(k) + dust%hflux_in(k) ) * dust_Fe_scavenge_scale + + (dust%sflux_in(k) + dust%hflux_in(k) + dust%prod(k)) * dust_Fe_scavenge_scale + + + if ((DOC_loc(k) .gt. c0) .and. (DOCr_loc(k) .gt. c0)) then + sinking_mass = sinking_mass * (((DOC_loc(k) * 2.0_r8) + DOCr_loc(k)) / DOCr_loc(k)) + endif + + Fe_scavenge_rate(k) = parm_Fe_scavenge_rate0 * sinking_mass Lig_scavenge_rate = parm_Lig_scavenge_rate0 * sinking_mass @@ -2572,7 +2649,7 @@ subroutine compute_large_detritus_prod(k, domain, marbl_tracer_indices, & type(marbl_tracer_index_type), intent(in) :: marbl_tracer_indices type(zooplankton_derived_terms_type), intent(in) :: zooplankton_derived_terms type(autotroph_derived_terms_type), intent(in) :: autotroph_derived_terms - real(r8), intent(in) :: Fe_scavenge + real(r8), intent(in) :: Fe_scavenge(domain%km) type(column_sinking_particle_type), intent(inout) :: POC type(column_sinking_particle_type), intent(inout) :: POP type(column_sinking_particle_type), intent(inout) :: PON @@ -2606,6 +2683,7 @@ subroutine compute_large_detritus_prod(k, domain, marbl_tracer_indices, & auto_loss_poc => autotroph_derived_terms%auto_loss_poc(:,k), & ! input remaining_P_pop => autotroph_derived_terms%remaining_P_pop(:,k), & ! input remaining_N_pon => autotroph_derived_terms%remaining_N_pon(:,k), & ! input + remaining_Fe_pfe => autotroph_derived_terms%remaining_Fe_pfe(:,k), & ! input zoo_loss_poc => zooplankton_derived_terms%zoo_loss_poc(:,k), & ! input zoo_graze_poc => zooplankton_derived_terms%zoo_graze_poc(:,k) & ! input ) @@ -2659,6 +2737,14 @@ subroutine compute_large_detritus_prod(k, domain, marbl_tracer_indices, & DON_loss_N_bal = c0 endif + !----------------------------------------------------------------------- + ! Iron + !----------------------------------------------------------------------- + + P_iron%prod(k) = (sum(zoo_loss_poc(:)) + sum(zoo_graze_poc(:))) * Qfe_zoo & + + Fe_scavenge(k) + sum(remaining_Fe_pfe(:)) + + !----------------------------------------------------------------------- ! large detrital CaCO3 ! 33% of CaCO3 is remin when phyto are grazed @@ -2667,7 +2753,8 @@ subroutine compute_large_detritus_prod(k, domain, marbl_tracer_indices, & P_CaCO3%prod(k) = c0 do auto_ind = 1, autotroph_cnt if (marbl_tracer_indices%auto_inds(auto_ind)%CaCO3_ind > 0) then - P_CaCO3%prod(k) = P_CaCO3%prod(k) + ((c1 - f_graze_CaCO3_remin) * auto_graze(auto_ind) & + P_CaCO3%prod(k) = P_CaCO3%prod(k) & + + ((c1 - f_graze_CaCO3_remin) * auto_graze(auto_ind) & + auto_loss(auto_ind) + auto_agg(auto_ind)) * QCaCO3(auto_ind) endif end do @@ -2688,21 +2775,11 @@ subroutine compute_large_detritus_prod(k, domain, marbl_tracer_indices, & end do !----------------------------------------------------------------------- - ! Dust + ! Dust prod now set in compute_scavenging !----------------------------------------------------------------------- - dust%prod(k) = c0 + ! dust%prod(k) = c0 - !----------------------------------------------------------------------- - ! Iron - !----------------------------------------------------------------------- - - P_iron%prod(k) = (sum(zoo_loss_poc(:)) + sum(zoo_graze_poc(:))) * Qfe_zoo + Fe_scavenge - - do auto_ind = 1, autotroph_cnt - P_iron%prod(k) = P_iron%prod(k) + Qfe(auto_ind) * & - (auto_agg(auto_ind) + auto_graze_poc(auto_ind) + auto_loss_poc(auto_ind)) - end do end associate end subroutine compute_large_detritus_prod @@ -2715,7 +2792,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & P_SiO2, dust, P_iron, QA_dust_def, & tracer_local, carbonate, denitrif_C_N, & sed_denitrif, other_remin, & - fesedflux, marbl_tracer_indices, & + fesedflux, feRedsedflux, feventflux, Lig_prod, marbl_tracer_indices, & glo_avg_fields_interior_tendency, marbl_status_log) ! Compute outgoing fluxes and remineralization terms. Assumes that @@ -2795,7 +2872,10 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & real (r8) , intent(in) :: p_remin_scalef ! Particulate Remin Scale Factor real (r8), dimension(:,:) , intent(in) :: tracer_local ! local copies of model tracer concentrations type(carbonate_type) , intent(in) :: carbonate - real(r8) , intent(in) :: fesedflux ! sedimentary Fe input + real (r8) , intent(in) :: fesedflux(domain%km) ! sedimentary Fe input + real (r8) , intent(in) :: feRedsedflux(domain%km) ! sedimentary Red Fe input + real (r8) , intent(in) :: feventflux(domain%km) ! vent Fe input + real (r8), dimension(:) , intent(inout) :: Lig_prod ! vent ligand inputs type(column_sinking_particle_type), intent(inout) :: POC ! base units = nmol C type(column_sinking_particle_type), intent(inout) :: POP ! base units = nmol P type(column_sinking_particle_type), intent(inout) :: PON ! base units = nmol N @@ -2807,7 +2887,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & real (r8) , intent(inout) :: QA_dust_def ! incoming deficit in the QA(dust) POC flux real (r8) , intent(inout) :: denitrif_C_N(:) ! C/N of denitrification real (r8), dimension(:) , intent(inout) :: sed_denitrif ! sedimentary denitrification (umolN/cm^2/s) - real (r8), dimension(:) , intent(inout) :: other_remin ! sedimentary remin not due to oxic or denitrification + real (r8), dimension(:) , intent(inout) :: other_remin ! sedimentary remin not oxic or denitrification type(marbl_particulate_share_type), intent(inout) :: marbl_particulate_share type(marbl_tracer_index_type) , intent(in) :: marbl_tracer_indices real (r8) , intent(inout) :: glo_avg_fields_interior_tendency(:) @@ -2841,8 +2921,8 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & flux_top, flux_bot, & ! total (sflux+hflux) particulate fluxes at top and bottom of layer wbot, & ! weight for interpolating fluxes vertically flux_alt, & ! flux to floor in alternative units, to match particular parameterizations - bury_frac, & ! fraction of flux hitting floor that gets buried - dz_loc, dzr_loc ! dz, dzr at a particular i, j location + bury_frac, & ! fraction of flux hitting floor that gets buried + dz_loc, dzr_loc, dzr_mod ! dz, dzr, dzr_mod at a particular i, j location real (r8) :: particulate_flux_ref_depth_cm @@ -2892,7 +2972,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & if (p_remin_scalef /= c1) scalelength = scalelength / p_remin_scalef DECAY_Hard = exp(-delta_z(k) * p_remin_scalef / 4.0e6_r8) - DECAY_HardDust = exp(-delta_z(k) * p_remin_scalef / 1.2e8_r8) + DECAY_HardDust = exp(-delta_z(k) * p_remin_scalef / 9.54641e7_r8) poc_error = .false. dz_loc = delta_z(k) @@ -2900,13 +2980,15 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & if (k <= column_kmt) then dzr_loc = c1 / dz_loc + dzr_mod = ((dz_loc * 0.01)**(-0.343)) + poc_diss = POC%diss sio2_diss = P_SiO2%diss caco3_diss = P_CaCO3%diss dust_diss = dust%diss !----------------------------------------------------------------------- - ! increase POC diss length scale where O2 concentrations are low + ! increase POC diss length scale where O2 concentrations are low !----------------------------------------------------------------------- if (O2_loc < o2_sf_o2_range_hi) then @@ -2941,7 +3023,6 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & ! The outoing fluxes for ballast materials are from the ! solution of the coresponding continuous ODE across the model ! level. The ODE has a constant source term and linear decay. - ! It is assumed that there is no sub-surface dust production. !----------------------------------------------------------------------- P_CaCO3%sflux_out(k) = P_CaCO3%sflux_in(k) * decay_CaCO3 + & @@ -2959,16 +3040,20 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & P_CaCO3_ALT_CO2%prod(k) * (P_CaCO3_ALT_CO2%gamma * dz_loc) P_SiO2%sflux_out(k) = P_SiO2%sflux_in(k) * decay_SiO2 + & - P_SiO2%prod(k) * ((c1 - P_SiO2%gamma) * (c1 - decay_SiO2) & - * sio2_diss) + P_SiO2%prod(k) * ((c1 - P_SiO2%gamma) * (c1 - decay_SiO2) * sio2_diss) P_SiO2%hflux_out(k) = P_SiO2%hflux_in(k) * DECAY_Hard + & P_SiO2%prod(k) * (P_SiO2%gamma * dz_loc) + + dust%sflux_out(k) = dust%sflux_in(k) * decay_dust + dust%hflux_in(k) = dust%hflux_in(k) + dust%prod(k) + dust%hflux_out(k) = dust%hflux_in(k) * DECAY_HardDust + !----------------------------------------------------------------------- ! Compute how much POC_PROD is available for deficit reduction ! and excess POC flux after subtracting off fraction of non-dust @@ -3040,8 +3125,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & POC_PROD_avail * ((c1 - decay_POC_E) * poc_diss) !----------------------------------------------------------------------- - ! Compute remineralization terms. It is assumed that there is no - ! sub-surface dust production. + ! Compute remineralization terms !----------------------------------------------------------------------- P_CaCO3%remin(k) = P_CaCO3%prod(k) + & @@ -3067,6 +3151,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & !----------------------------------------------------------------------- ! Compute iron remineralization and flux out. !----------------------------------------------------------------------- + P_iron%remin(k) = c0 if (POC%sflux_in(k) + POC%hflux_in(k) == c0) then P_iron%remin(k) = (POC%remin(k) * parm_Red_Fe_C) @@ -3076,10 +3161,6 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & (POC%sflux_in(k) + POC%hflux_in(k))) endif - ! add term for desorption of iron from sinking particles - P_iron%remin(k) = P_iron%remin(k) + & - (P_iron%sflux_in(k) * parm_Fe_desorption_rate0) - P_iron%sflux_out(k) = P_iron%sflux_in(k) + dz_loc * & ((c1 - P_iron%gamma) * P_iron%prod(k) - P_iron%remin(k)) @@ -3090,17 +3171,36 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & endif !----------------------------------------------------------------------- - ! Compute iron release from dust remin/dissolution and add in Fe source - ! from sediments + ! Compute iron release from dust remin/dissolution and other Fe sources !----------------------------------------------------------------------- + !----------------------------------------------------------------------- + ! add term for desorption of iron from sinking particles + ! dzr_mod increases slowly relative to dzr_loc + ! it accounts for the increasing sinking speed of particles with depth + ! less desorption with depth (as mean sinking speed increases) + !----------------------------------------------------------------------- + + P_iron%remin(k) = P_iron%remin(k) + & + (P_iron%sflux_in(k) * parm_Fe_desorption_rate0 * dzr_mod) + + + P_iron%remin(k) = P_iron%remin(k) & + + (dust%remin(k) * dust_to_Fe) & + + (fesedflux(k) * dzr_loc) & + + (feRedsedflux(k) * dzr_loc) & + + (feventflux(k) * dzr_loc) - P_iron%remin(k) = P_iron%remin(k) & - + dust%remin(k) * dust_to_Fe & - + (fesedflux * dzr_loc) P_iron%hflux_out(k) = P_iron%hflux_in(k) + + !----------------------------------------------------------------------- + ! add ligand source from vents, proportional to Fe input + !----------------------------------------------------------------------- + + Lig_prod(k) = feventflux(k) * dzr_loc * 0.22 + !------------------------------------------------------------------------ ! compute POP and PON remin and flux out, following code for iron !------------------------------------------------------------------------ @@ -3294,7 +3394,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & endif sed_denitrif(1:k) = bot_flux_to_tend(1:k) * parm_sed_denitrif_coeff * POC%to_floor & - * (0.06_r8 + 0.19_r8 * 0.99_r8**(O2_loc-NO3_loc)) + * (0.06_r8 + 0.19_r8 * 0.99_r8**(O2_loc-NO3_loc)) * 1.5_r8 flux_alt = POC%to_floor*1.0e-6_r8*spd*365.0_r8 ! convert to mmol/cm^2/year other_remin(1:k) = min(bot_flux_to_tend(1:k) * & @@ -3409,7 +3509,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, & endif !----------------------------------------------------------------------- ! Remove all Piron and dust that hits bottom, sedimentary Fe source - ! accounted for by fesedflux elsewhere. Remin C, N, P not buried. + ! accounted for elsewhere. Remin C, N, P not buried. !----------------------------------------------------------------------- P_iron%to_floor = P_iron%sflux_out(k) + P_iron%hflux_out(k) @@ -3472,10 +3572,10 @@ subroutine compute_Lig_terms(km, PAR_nsubcols, marbl_tracer_indices, & ! and production of DOM !----------------------------------------------------------------------- - Lig_prod(:) = (POC_remin(:) * remin_to_Lig) + (DOC_prod(:) * remin_to_Lig) + Lig_prod(:) = Lig_prod(:) + (POC_remin(:) * remin_to_Lig) + (DOC_prod(:) * remin_to_Lig) !---------------------------------------------------------------------- - ! Ligand losses due to photochemistry in first ocean layer + ! Ligand losses due to photochemistry in first ocean layer 1/3months ! ligand photo-oxidation a function of PAR (2/5/2015) !---------------------------------------------------------------------- @@ -3485,14 +3585,14 @@ subroutine compute_Lig_terms(km, PAR_nsubcols, marbl_tracer_indices, & rate_per_sec = c0 do subcol_ind = 1, PAR_nsubcols if ((PAR_col_frac(subcol_ind) > c0) .and. (PAR%interface(0,subcol_ind) > 1.0_r8)) then - rate_per_sec_subcol = (log(PAR%interface(0,subcol_ind))*0.4373_r8)*(10.0e2_r8/dz1)*((12.0_r8/3.0_r8)*yps) ! 1/(3 months) + rate_per_sec_subcol = (log(PAR%interface(0,subcol_ind))*0.4373_r8)*(10.0e2_r8/dz1)*(0.009_r8*dps) rate_per_sec = rate_per_sec + PAR_col_frac(subcol_ind) * rate_per_sec_subcol endif end do Lig_photochem(1) = Lig_loc(1) * rate_per_sec !---------------------------------------------------------------------- - ! Ligand losses due to uptake/degradation (by hetrotrophic bacteria) + ! Ligand losses due to uptake/degradation (by hetrotrophic bacteria), scaled by POCremin !---------------------------------------------------------------------- Lig_deg(:) = POC_remin(:) * parm_Lig_degrade_rate0 @@ -3501,7 +3601,7 @@ subroutine compute_Lig_terms(km, PAR_nsubcols, marbl_tracer_indices, & ! total ligand loss !---------------------------------------------------------------------- - Lig_loss(:) = Lig_scavenge(:) + 0.20_r8*sum(photoFe(:,:),dim=1) + Lig_photochem(:) + Lig_deg(:) + Lig_loss(:) = Lig_scavenge(:) + 0.25_r8*sum(photoFe(:,:),dim=1) + Lig_photochem(:) + Lig_deg(:) end associate @@ -3603,11 +3703,12 @@ subroutine compute_denitrif(km, marbl_tracer_indices, tracer_local, DOC_remin, D - other_remin(k)) / denitrif_C_N(k) - sed_denitrif(k)) ! scale down denitrif if computed rate would consume all NO3 in 10 days - if (NO3_loc(k) < ((c10*spd)*(denitrif(k)+sed_denitrif(k)))) then - work = NO3_loc(k) / ((c10*spd)*(denitrif(k)+sed_denitrif(k))) - denitrif(k) = denitrif(k) * work - sed_denitrif(k) = sed_denitrif(k) * work - end if + ! if (NO3_loc(k) < ((c10*spd)*(denitrif(k)+sed_denitrif(k)))) then + ! work = NO3_loc(k) / ((c10*spd)*(denitrif(k)+sed_denitrif(k))) + ! denitrif(k) = denitrif(k) * work + ! sed_denitrif(k) = sed_denitrif(k) * work + ! end if + end do end associate @@ -3654,6 +3755,8 @@ subroutine compute_local_tendencies(km, marbl_tracer_indices, autotroph_derived_ integer :: k, auto_ind, zoo_ind, n real(r8) :: auto_sum real(r8) :: o2_production_denom + + !----------------------------------------------------------------------- associate( & @@ -3682,8 +3785,10 @@ subroutine compute_local_tendencies(km, marbl_tracer_indices, autotroph_derived_ CaCO3_form => autotroph_derived_terms%CaCO3_form(:,:), & ! prod. of CaCO3 by small phyto (mmol CaCO3/m^3/sec) Nfix => autotroph_derived_terms%Nfix(:,:), & ! total Nitrogen fixation (mmol N/m^3/sec) Nexcrete => autotroph_derived_terms%Nexcrete(:,:), & ! fixed N excretion by diazotrophs - remaining_P_dip => autotroph_derived_terms%remaining_P_dip(:,:), & ! remaining_P from mort routed to remin - remaining_N_din => autotroph_derived_terms%remaining_N_din(:,:), & ! remaining_N from mort routed to remin + remaining_P_dip => autotroph_derived_terms%remaining_P_dip(:,:), & ! remaining_P from grazing mort routed to remin + remaining_N_din => autotroph_derived_terms%remaining_N_din(:,:), & ! remaining_N from grazing mort routed to remin + remaining_Fe_dfe => autotroph_derived_terms%remaining_Fe_dfe(:,:), & ! remaining_Fe from grazing mort routed to remin + remaining_Fe_pfe => autotroph_derived_terms%remaining_Fe_pfe(:,:), & ! remaining_Fe from phyto mort routed to particulate sinking x_graze_zoo => zooplankton_derived_terms%x_graze_zoo(:,:), & ! {auto, zoo}_graze routed to zoo (mmol C/m^3/sec) zoo_graze => zooplankton_derived_terms%zoo_graze(:,:), & ! zooplankton losses due to grazing (mmol C/m^3/sec) @@ -3733,11 +3838,6 @@ subroutine compute_local_tendencies(km, marbl_tracer_indices, autotroph_derived_ interior_tendencies(no3_ind,k) = nitrif(k) - denitrif(k) - sed_denitrif(k) - sum(NO3_V(:,k)) - !interior_tendencies(nh4_ind,k) = -sum(NH4_V(:,k)) - nitrif(k) + DON_remin(k) + DONr_remin(k) & - ! + Qn_zoo * (sum(zoo_loss_dic(:,k)) + sum(zoo_graze_dic(:,k)) + sum(auto_loss_dic(:,k)) & - ! + sum(auto_graze_dic(:,k)) + DOC_prod(k)*(c1 - f_toDON)) & - ! + PON_remin(k) * (c1 - PONremin_refract) - interior_tendencies(nh4_ind,k) = DON_remin(k) + DONr_remin(k) - sum(NH4_V(:,k)) - nitrif(k) & + ((c1 - PONremin_refract) * PON_remin(k)) + sum(remaining_N_din(:,k)) & @@ -3754,23 +3854,17 @@ subroutine compute_local_tendencies(km, marbl_tracer_indices, autotroph_derived_ !----------------------------------------------------------------------- interior_tendencies(po4_ind,k) = DOP_remin(k) + DOPr_remin(k) - sum(PO4_V(:,k)) & - + (c1 - POPremin_refract) * POP_remin(k) + sum(remaining_P_dip(:,k)) & - + Qp_zoo * (sum(zoo_loss_dic(:,k)) + sum(zoo_graze_dic(:,k))) + + (c1 - POPremin_refract) * POP_remin(k) + sum(remaining_P_dip(:,k)) & + + Qp_zoo * (sum(zoo_loss_dic(:,k)) + sum(zoo_graze_dic(:,k))) !----------------------------------------------------------------------- ! dissolved iron !----------------------------------------------------------------------- interior_tendencies(fe_ind,k) = P_iron_remin(k) - sum(photoFe(:,k)) - Fe_scavenge(k) & - + Qfe_zoo * (sum(zoo_loss_dic(:,k)) + sum(zoo_loss_doc(:,k)) & - + sum(zoo_graze_dic(:,k)) + sum(zoo_graze_doc(:,k))) - - do auto_ind = 1, autotroph_cnt - interior_tendencies(fe_ind,k) = interior_tendencies(fe_ind,k) & - + (Qfe(auto_ind,k) * (auto_loss_dic(auto_ind,k) + auto_graze_dic(auto_ind,k))) & - + auto_graze_zootot(auto_ind,k) * (Qfe(auto_ind,k) - Qfe_zoo) & - + (Qfe(auto_ind,k) * (auto_loss_doc(auto_ind,k) + auto_graze_doc(auto_ind,k))) - end do + + sum(remaining_Fe_dfe(:,k)) & + + Qfe_zoo * (sum(zoo_loss_dic(:,k)) + sum(zoo_loss_doc(:,k)) & + + sum(zoo_graze_dic(:,k)) + sum(zoo_graze_doc(:,k))) !----------------------------------------------------------------------- ! iron binding ligand @@ -3787,9 +3881,8 @@ subroutine compute_local_tendencies(km, marbl_tracer_indices, autotroph_derived_ do auto_ind = 1, autotroph_cnt if (marbl_tracer_indices%auto_inds(auto_ind)%Si_ind > 0) then interior_tendencies(sio3_ind,k) = interior_tendencies(sio3_ind,k) - photoSi(auto_ind,k) & - + Qsi(auto_ind,k) * (f_graze_si_remin * auto_graze(auto_ind,k) & - + (c1 - autotroph_settings(auto_ind)%loss_poc) & - * auto_loss(auto_ind,k)) + + Qsi(auto_ind,k) * (f_graze_si_remin * auto_graze(auto_ind,k) & + + (c1 - autotroph_settings(auto_ind)%loss_poc) * auto_loss(auto_ind,k)) end if end do diff --git a/src/marbl_settings_mod.F90 b/src/marbl_settings_mod.F90 index 0c707196..7e5e8e44 100644 --- a/src/marbl_settings_mod.F90 +++ b/src/marbl_settings_mod.F90 @@ -96,19 +96,20 @@ module marbl_settings_mod real(r8), parameter :: & Q_10 = 1.7_r8, & ! factor for temperature dependence (non-dim) xkw_coeff = 6.97e-9_r8, & ! in s/cm, from a = 0.251 cm/hr s^2/m^2 in Wannikhof 2014 - parm_Red_D_C_P = 117.0_r8, & ! carbon:phosphorus - parm_Red_D_N_P = 16.0_r8, & ! nitrogen:phosphorus - parm_Red_D_O2_P = 170.0_r8, & ! oxygen:phosphorus - parm_Remin_D_O2_P = 138.0_r8, & ! oxygen:phosphorus - parm_Red_P_C_P = parm_Red_D_C_P, & ! carbon:phosphorus - parm_Red_D_C_N = parm_Red_D_C_P/parm_Red_D_N_P, & ! carbon:nitrogen - parm_Red_P_C_N = parm_Red_D_C_N, & ! carbon:nitrogen - parm_Red_D_C_O2 = parm_Red_D_C_P/parm_Red_D_O2_P, & ! carbon:oxygen - parm_Remin_D_C_O2 = parm_Red_D_C_P/parm_Remin_D_O2_P, & ! carbon:oxygen - parm_Red_P_C_O2 = parm_Red_D_C_O2, & ! carbon:oxygen + + parm_Red_D_C_P = 112.0_r8, & ! carbon:phosphorus + parm_Red_D_N_P = 16.0_r8, & ! nitrogen:phosphorus + parm_Remin_D_O2_P = 112.0_r8, & ! oxygen:phosphorus + parm_Red_D_O2_P = (parm_Remin_D_O2_P + 32.0_r8), & ! oxygen:phosphorus + parm_Red_P_C_P = parm_Red_D_C_P, & ! carbon:phosphorus + parm_Red_D_C_N = parm_Red_D_C_P/parm_Red_D_N_P, & ! carbon:nitrogen + parm_Red_P_C_N = parm_Red_D_C_N, & ! carbon:nitrogen + parm_Red_D_C_O2 = parm_Red_D_C_P/parm_Red_D_O2_P, & ! carbon:oxygen ** + parm_Remin_D_C_O2 = parm_Red_D_C_P/parm_Remin_D_O2_P, & ! carbon:oxygen ** + parm_Red_P_C_O2 = parm_Red_D_C_O2, & ! carbon:oxygen ** parm_Red_Fe_C = 3.0e-6_r8, & ! iron:carbon - parm_Red_D_C_O2_diaz = parm_Red_D_C_P/150.0_r8 ! carbon:oxygen - ! for diazotrophs + parm_Red_D_C_O2_diaz = parm_Red_D_C_P/(parm_Remin_D_O2_P + 2.0_r8)! carbon:oxygen ** + ! for diazotrophs ! Misc. Rate constants real(r8), parameter :: & @@ -130,11 +131,12 @@ module marbl_settings_mod caco3_poc_min = 0.40_r8, & ! minimum proportionality between ! QCaCO3 and grazing losses to POC ! (mmol C/mmol CaCO3) - spc_poc_fac = 0.13_r8, & ! small phyto grazing factor (1/mmolC) - f_graze_sp_poc_lim = 0.36_r8, & + spc_poc_fac = 0.04_r8, & ! small phyto grazing factor (1/mmolC) + f_graze_sp_poc_lim = 0.33_r8, & f_photosp_CaCO3 = 0.40_r8, & ! proportionality between small phyto + ! production and CaCO3 production - f_graze_si_remin = 0.50_r8, & ! fraction of diatom Si grazing which is remin + f_graze_si_remin = 0.22_r8, & ! fraction of diatom Si grazing which is remin f_toDON = 0.05_r8, & ! fraction of reamining_N to DON f_toDOP = 0.15_r8 ! fraction of remaining_P to DOP @@ -143,8 +145,9 @@ module marbl_settings_mod ! SET parmaeters and RATIOS for N/C, P/C, SiO3/C, Fe/C, etc... real(r8), parameter :: & - Q = 16.0_r8 / 117.0_r8, & ! N/C ratio (mmol/mmol) of phyto & zoo - Qfe_zoo = 3.0e-6_r8 ! zooplankton Fe/C ratio + Q = 16.0_r8 / 112.0_r8, & ! N/C ratio (mmol/mmol) of phyto & zoo + Qfe_zoo = 6.0e-6_r8 ! zooplankton Fe/C ratio + ! parameters in GalbraithMartiny Pquota Model^M !POpt = 0.6_r8, & !gQp_max = 0.01111111_r8, & ! maximum NP = 16 @@ -154,10 +157,12 @@ module marbl_settings_mod ! loss term threshold parameters, chl:c ratios real(r8), parameter :: & - thres_z1_auto = 80.0e2_r8, & ! autotroph threshold = C_loss_thres for z shallower than this (cm) - thres_z2_auto = 120.0e2_r8, & ! autotroph threshold = 0 for z deeper than this (cm) - thres_z1_zoo = 110.0e2_r8, & ! zooplankton threshold = C_loss_thres for z shallower than this (cm) - thres_z2_zoo = 150.0e2_r8, & ! zooplankton threshold = 0 for z deeper than this (cm) + + thres_z1_auto = 50.0e2_r8, & ! autotroph threshold = C_loss_thres for z shallower than this (cm) + thres_z2_auto = 100.0e2_r8, & ! autotroph threshold = 0 for z deeper than this (cm) + thres_z1_zoo = 40.0e2_r8, & ! zooplankton threshold = C_loss_thres for z shallower than this (cm) + thres_z2_zoo = 80.0e2_r8, & ! zooplankton threshold = 0 for z deeper than this (cm) + CaCO3_temp_thres1 = 4.0_r8, & ! upper temp threshold for CaCO3 prod CaCO3_temp_thres2 = -2.0_r8, & ! lower temp threshold CaCO3_sp_thres = 2.5_r8 ! bloom condition thres (mmolC/m3) @@ -3035,3 +3040,4 @@ end subroutine set_derived_from_temp_func_form !***************************************************************************** end module marbl_settings_mod +