Changes to pass call_compute_subroutines tests

Modify marbl_io_mod to initialize N tracers to C tracers * 16/112 and set
feredsedflux and feventflux = 0.
Modify marbl_diagnostics_mod.F90 to include dust%remin*dust_to_Fe in integrand,
not righthand-side of iron conservation check
Modify marbl_interior_tendency_mod.F90 to sflux_in * desportion in P_iron%remin
before computing sflux_out.

Also cleaned up a few small things in interior_tendency_mod (use unit_system
instead of assuming cgs; initialize Lig_prod and declare it as inout in some
routines)

src/marbl_diagnostics_mod.F90

@@ -4681,20 +4681,21 @@ subroutine store_diagnostics_iron_fluxes(marbl_domain, P_iron, dust, &
     ! 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))
+           + (Qfe_zoo * sum(interior_tendencies(marbl_tracer_indices%zoo_inds(:)%C_ind, :),dim=1)) &
+           - (dust%remin * dust_to_Fe)
 
 
     call marbl_diagnostics_share_compute_vertical_integrals(work, delta_z, kmt, unit_system, &
       full_depth_integral=diags(ind%Jint_Fetot)%field_2d(1),                     &
-      integrated_terms = P_iron%sed_loss - fesedflux - feRedsedflux - feventflux - (dust%remin * dust_to_Fe))
+      integrated_terms = P_iron%sed_loss - fesedflux - feRedsedflux - feventflux)
 
 
    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
+      return
    end if
 
     end associate

src/marbl_interior_tendency_mod.F90

@@ -1316,26 +1316,22 @@ subroutine compute_autotroph_elemental_ratios(km, autotroph_local, marbl_tracer_
         !------------------------------------------------------------------------
         if (lvariable_PtoC) then
           Qp(auto_ind,:) = auto_P(auto_ind,:) / (auto_C(auto_ind,:) + epsC)
-
-          !!-----------------------------------------------------------------------
-          !!-- Calculate Qp for new growth based on Galbraith and Martiny (2015), with min. N/P
-          !! - 14= 0.00976801, 14.5 = 0.00944239 15= 0.00911677 15.5=0.00882272 16= 0.00854701
-          !! - std intercept 6.0 = 166.66maxCP, 5.26=190, 4.0 = 250, 3.0 = 333.33 (commented out)
-          !!-----------------------------------------------------------------------
-          !gQp(auto_ind,:) = min((((PquotaSlope * PO4_loc(:)) + PquotaIntercept) * 0.001_r8), PquotaMinNP)
-
           gQp(auto_ind,:) = autotroph_settings(auto_ind)%gQp_max
 
           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,:) * WORK1 / autotroph_settings(auto_ind)%NOpt, &
-                autotroph_settings(auto_ind)%gQp_min * 0.55_r8)
-          endwhere
+
+          ! Apply threshold if using variable N:C
+          if (lvariable_NtoC) then
+            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.55_r8)
+            endwhere
+          endif
 
         else
           Qp(auto_ind,:) = autotroph_settings(auto_ind)%Qp_fixed
@@ -3024,7 +3020,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, p_remin_scalef
      if (k <= column_kmt) then
 
         dzr_loc    = c1 / dz_loc
-        dzr_mod = ((dz_loc * 0.01_r8)**(-0.343_r8))
+        dzr_mod = ((dz_loc * unit_system%len2m)**(-0.343_r8))
 
         poc_diss   = POC%diss
         sio2_diss  = P_SiO2%diss
@@ -3205,6 +3201,16 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, p_remin_scalef
                 (POC%sflux_in(k) + POC%hflux_in(k)))
         endif
 
+        !-----------------------------------------------------------------------
+        ! 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%sflux_out(k) = P_iron%sflux_in(k) + dz_loc * &
              ((c1 - P_iron%gamma) * P_iron%prod(k) - P_iron%remin(k))
 
@@ -3218,17 +3224,6 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, p_remin_scalef
         !  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)        &
@@ -3399,6 +3394,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, p_remin_scalef
      P_iron%sed_loss(k)          = c0
      dust%sed_loss(k)            = c0
      denitrif_C_N(k)             = 0.9375_r8
+     Lig_prod(k)                 = c0
 
      if (k == column_kmt) then
 
@@ -3594,7 +3590,7 @@ subroutine compute_Lig_terms(km, PAR_nsubcols, marbl_tracer_indices, &
     real(r8),                      intent(in)  :: Lig_scavenge(km)
     real(r8),                      intent(in)  :: photoFe(autotroph_cnt,km)
     type(unit_system_type),        intent(in)  :: unit_system
-    real(r8),                      intent(out) :: Lig_prod(km)
+    real(r8),                      intent(inout) :: Lig_prod(km)
     real(r8),                      intent(out) :: Lig_photochem(km)
     real(r8),                      intent(out) :: Lig_deg(km)
     real(r8),                      intent(out) :: Lig_loss(km)

tests/driver_src/marbl_io_mod.F90

@@ -485,6 +485,12 @@ subroutine marbl_io_read_forcing_field(col_id, lat, unit_system_opt, forcing_fie
           forcing_fields_out(n)%field_0d(col_id) = -5.8_r8
         end if
         cycle
+      else if (trim(varname) == 'iron_redsed_flux') then
+        forcing_fields_out(n)%field_1d(col_id,:) = 0._r8
+        cycle
+      else if (trim(varname) == 'iron_vent_flux') then
+        forcing_fields_out(n)%field_1d(col_id,:) = 0._r8
+        cycle
       end if
 
       ! Get netcdf varid
@@ -552,6 +558,18 @@ subroutine get_tracer_name_file_and_scale_factor(tracer_name, tracer_name_file,
       return
     end if
     select case (trim(tracer_name))
+      case('spN')
+        tracer_name_file = 'spC'
+        scale_factor = 16._r8 / 112._r8
+      case('diatN')
+        tracer_name_file = 'diatC'
+        scale_factor = 16._r8 / 112._r8
+      case('diazN')
+        tracer_name_file = 'diazC'
+        scale_factor = 16._r8 / 112._r8
+      case('coccoN')
+        tracer_name_file = 'coccoC'
+        scale_factor = 16._r8 / 112._r8
       case('coccoChl')
         tracer_name_file = 'spChl'
         scale_factor = 0.07_r8
@@ -1318,6 +1336,18 @@ subroutine get_forcing_varname_rank_and_conv_factor(forcing_name, unit_system_op
         if (trim(unit_system_opt) == 'mks') &
           ! convert from nmol/cm^2/s -> mmol/m^2/s
           conv_factor = 0.01_r8
+      case('Iron Red Sediment Flux')
+        varname = 'iron_redsed_flux'
+        rank = 1
+        if (trim(unit_system_opt) == 'mks') &
+          ! convert from nmol/cm^2/s -> mmol/m^2/s
+          conv_factor = 0.01_r8
+      case('Iron Vent Flux')
+        varname = 'iron_vent_flux'
+        rank = 1
+        if (trim(unit_system_opt) == 'mks') &
+          ! convert from nmol/cm^2/s -> mmol/m^2/s
+          conv_factor = 0.01_r8
       case('O2 Consumption Scale Factor')
         varname = 'o2_consumption_scalef'
         rank = 1