First commit following code review

Keith and I reviewed the fortran changes in src/ and this commit addresses
concerns raised during that review. Most of it is straight forward, I do want
to highlight one change that I think still needs improvement:

Instead of passing unit_system_opt ('cgs' or 'mks') to
marbl_instance%surface_flux_output%add_output(), we call
marbl_instance%get_conc_flux_units and then pass the conc_flux_units directly.
This still feels kludgy, but the way the interface and the
marbl_output_for_GCM_type are defined it was the only option for a quick fix.
I'd like to revisit this interface and come up with a complete redesign of the
way we register fields the GCM would like MARBL to return.

MARBL_tools/netcdf_comparison.py

@@ -150,13 +150,14 @@ def _reduce_to_matching_variables(ds_base, ds_new):
 def _get_conversion_factor(ds_base, ds_new, var):
  unit_conversion = {key: {} for key in
  ['cm/s', 'cm','nmol/cm^3', 'nmol/cm^3/s', 'nmol/cm^2/s',
- 'g/cm^3/s', 'g/cm^2/s', 'meq/m^3', 'meq/m^3/s',
+ '(nmol/cm^3)^-1 s^-1', 'g/cm^3/s', 'g/cm^2/s', 'meq/m^3', 'meq/m^3/s',
  'neq/cm^2/s', 'meq/m^3 cm/s', 'mg/m^3 cm/s']}
  unit_conversion['cm/s']['m/s'] = 0.01 # cm/s -> m/s
  unit_conversion['cm']['m'] = 0.01 # cm -> m
  unit_conversion['nmol/cm^3']['mmol/m^3'] = 1 # nmol/cm^3 -> mmol/m^3
  unit_conversion['nmol/cm^3/s']['mmol/m^3/s'] = 1 # nmol/cm^3/s -> mmol/m^3/s
  unit_conversion['nmol/cm^2/s']['mmol/m^2/s'] = 0.01 # nmol/cm^2/s -> mmol/m^2/s
+ unit_conversion['(nmol/cm^3)^-1 s^-1']['(mmol/m^3)^-1 s^-1'] = 1 # same as nmol/cm^3 -> mmol/m^3
  unit_conversion['g/cm^3/s']['kg/m^3/s'] = 1000 # g/cm^3/s -> kg/m^3/s
  unit_conversion['g/cm^2/s']['kg/m^2/s'] = 10 # g/cm^2/s -> kg/m^2/s
  unit_conversion['meq/m^3']['neq/cm^3'] = 1 # meq/m^3 -> neq/cm^3

src/marbl_ciso_diagnostics_mod.F90

@@ -793,7 +793,7 @@ subroutine marbl_ciso_diagnostics_init( &
  ! TODO: generalize units for cgs or mks
  lname = trim(autotroph_settings(n)%lname) // ' instanteous growth rate over [CO2*]'
  sname = 'CISO_mui_to_co2star_' // trim(autotroph_settings(n)%sname)
- units = 'm^3/mmol/s'
+ write(units, '(3A)') '(', trim(unit_system%conc_units), ')^-1 s^-1'
  vgrid = 'layer_avg'
  truncate = .false.
  call diags%add_diagnostic(lname, sname, units, vgrid, truncate, &

src/marbl_ciso_init_mod.F90

@@ -14,19 +14,19 @@ module marbl_ciso_init_mod
 
  !*****************************************************************************
 
- subroutine marbl_ciso_init_tracer_metadata(unit_system,  &
- marbl_tracer_metadata, &
- marbl_tracer_indices)
+ subroutine marbl_ciso_init_tracer_metadata(unit_system, &
+ marbl_tracer_indices, &
+ marbl_tracer_metadata)
 
  ! Set tracer and forcing metadata
  use marbl_settings_mod, only : lecovars_full_depth_tavg
  use marbl_settings_mod, only : autotroph_cnt
  use marbl_settings_mod, only : autotroph_settings
  use marbl_settings_mod, only : unit_system_type
 
- type(unit_system_type),  intent(in) :: unit_system
- type (marbl_tracer_metadata_type) , intent(inout) :: marbl_tracer_metadata(:) ! descriptors for each tracer
- type(marbl_tracer_index_type) , intent(in) :: marbl_tracer_indices
+ type(unit_system_type), intent(in) :: unit_system
+ type(marbl_tracer_index_type), intent(in)  :: marbl_tracer_indices
+ type(marbl_tracer_metadata_type), intent(inout) :: marbl_tracer_metadata(:) ! descriptors for each tracer
 
  !-----------------------------------------------------------------------
  ! local variables

src/marbl_ciso_interior_tendency_mod.F90

@@ -1077,6 +1077,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, tracer_local,
  !----------------------------------------------------------------------------------------
 
  use marbl_constants_mod, only : spd
+ use marbl_constants_mod, only : spy
  use marbl_settings_mod , only : denitrif_C_N
  use marbl_settings_mod , only : parm_sed_denitrif_coeff
  use marbl_settings_mod , only : caco3_bury_thres_iopt
@@ -1261,7 +1262,7 @@ subroutine compute_particulate_terms(k, domain, bot_flux_to_tend, tracer_local,
  sed_denitrif(1:k) = bot_flux_to_tend(1:k) * parm_sed_denitrif_coeff * POC_ciso%to_floor * &
  (0.06_r8 + 0.19_r8 * 0.99_r8**(O2_loc-NO3_loc))
 
- flux_alt = POC_ciso%to_floor*(unit_system%len2cm * (mpercm**3))*spd*365.0_r8 ! convert to mmol/cm^2/year
+ flux_alt = POC_ciso%to_floor*(unit_system%conc_flux2mmol_m2s * (mpercm**2))*spy ! convert to mmol/cm^2/year
  other_remin(1:k) = min(bot_flux_to_tend(1:k) * &
  min(0.1_r8 + flux_alt,0.5_r8) * (POC_ciso%to_floor - POC_ciso%sed_loss(k)), &
  bot_flux_to_tend(1:k) * (POC_ciso%to_floor - POC_ciso%sed_loss(k)) - &

src/marbl_co2calc_mod.F90

@@ -132,13 +132,13 @@ subroutine marbl_co2calc_surface( &
  )
 
  !---------------------------------------------------------------------------
- ! set unit conversion factors
+ ! set unit conversion factors (literature assumes conc in mol/kg)
  ! mass_to_vol: mol / kg -> mmol / m^3
  ! => mmol / mol * kg / m^3
  ! rho_sw is M / L^3 => rho_sw * 1e3 * kg/M * (L/m)^3
  !---------------------------------------------------------------------------
 
- mass_to_vol = rho_sw * (1e3 * unit_system%mass2kg * (unit_system%m2len)**3)
+ mass_to_vol = rho_sw * (1.e3_r8 * unit_system%mass2kg * (unit_system%m2len)**3)
  vol_to_mass = c1 / mass_to_vol
 
  !---------------------------------------------------------------------------
@@ -277,13 +277,13 @@ subroutine marbl_co2calc_interior(&
  )
 
  !---------------------------------------------------------------------------
- ! set unit conversion factors
+ ! set unit conversion factors (literature assumes conc in mol/kg)
  ! mass_to_vol: mol / kg -> mmol / m^3
  ! => mmol / mol * kg / m^3
  ! rho_sw is M / L^3 => rho_sw * 1e3 * kg/M * (L/m)^3
  !---------------------------------------------------------------------------
 
- mass_to_vol = rho_sw * (1e3 * unit_system%mass2kg * (unit_system%m2len)**3)
+ mass_to_vol = rho_sw * (1.e3_r8 * unit_system%mass2kg * (unit_system%m2len)**3)
  vol_to_mass = c1 / mass_to_vol
 
  !------------------------------------------------------------------------
@@ -708,13 +708,13 @@ subroutine comp_htotal(num_elements, num_active_elements, dic_in, ta_in, pt_in,
  )
 
  !---------------------------------------------------------------------------
- ! set unit conversion factors
+ ! set unit conversion factors (literature assumes conc in mol/kg)
  ! mass_to_vol: mol / kg -> mmol / m^3
  ! => mmol / mol * kg / m^3
  ! rho_sw is M / L^3 => rho_sw * 1e3 * kg/M * (L/m)^3
  !---------------------------------------------------------------------------
 
- mass_to_vol = rho_sw * (1e3 * unit_system%mass2kg * (unit_system%m2len)**3)
+ mass_to_vol = rho_sw * (1.e3_r8 * unit_system%mass2kg * (unit_system%m2len)**3)
  vol_to_mass = c1 / mass_to_vol
 
  !---------------------------------------------------------------------------
@@ -1126,13 +1126,13 @@ subroutine marbl_co2calc_co3_sat_vals(&
  !---------------------------------------------------------------------------
 
  !---------------------------------------------------------------------------
- ! set unit conversion factors
+ ! set unit conversion factors (literature assumes conc in mol/kg)
  ! mass_to_vol: mol / kg -> mmol / m^3
  ! => mmol / mol * kg / m^3
  ! rho_sw is M / L^3 => rho_sw * 1e3 * kg/M * (L/m)^3
  !---------------------------------------------------------------------------
 
- mass_to_vol = rho_sw * (1e3 * unit_system%mass2kg * (unit_system%m2len)**3)
+ mass_to_vol = rho_sw * (1.e3_r8 * unit_system%mass2kg * (unit_system%m2len)**3)
 
  salt_lim = max(salt(:),salt_min)
  tk = T0_Kelvin + temp(:)

src/marbl_diagnostics_mod.F90

@@ -1708,7 +1708,6 @@ subroutine marbl_diagnostics_init( &
 
  lname = 'Nitrification'
  sname = 'NITRIF'
- write(units, "(2A)") trim(unit_system%conc_units), '/s'
  units = unit_system%conc_tend_units
  vgrid = 'layer_avg'
  truncate = .false.
@@ -3136,7 +3135,7 @@ subroutine marbl_diagnostics_interior_tendency_compute( &
  real (r8), intent(in) :: column_o2(:)
  real (r8), intent(in) :: o2_production(:)
  real (r8), intent(in) :: o2_consumption(:)
- real (r8), intent(in) :: fe_scavenge_rate(domain%km) ! annual scavenging rate of iron as % of ambient
+ real (r8), intent(in) :: fe_scavenge_rate(domain%km) ! scavenging rate of iron as % of ambient
  real (r8), intent(in) :: fe_scavenge(domain%km) ! loss of dissolved iron, scavenging (mmol Fe/m^3/sec)
  real (r8), intent(in) :: Lig_prod(domain%km)
  real (r8), intent(in) :: Lig_loss(domain%km)

src/marbl_init_mod.F90

@@ -181,8 +181,8 @@ subroutine marbl_init_tracers(num_levels, &
  allocate(tracer_restore_vars(tracer_indices%total_cnt))
 
  ! Set up tracer metadata
- call marbl_init_tracer_metadata(tracer_metadata, tracer_indices, unit_system)
- call marbl_ciso_init_tracer_metadata(unit_system, tracer_metadata, tracer_indices)
+ call marbl_init_tracer_metadata(unit_system, tracer_indices, tracer_metadata)
+ call marbl_ciso_init_tracer_metadata(unit_system, tracer_indices, tracer_metadata)
 
  ! Log what tracers are being used
  call marbl_status_log%log_header('MARBL Tracer indices', subname)
@@ -211,15 +211,15 @@ end subroutine marbl_init_tracers
 
  !***********************************************************************
 
- subroutine marbl_init_tracer_metadata(marbl_tracer_metadata, marbl_tracer_indices, unit_system)
+ subroutine marbl_init_tracer_metadata(unit_system, marbl_tracer_indices, marbl_tracer_metadata)
 
  ! Set tracer and forcing metadata
 
  use marbl_settings_mod, only : lecovars_full_depth_tavg
 
- type (marbl_tracer_metadata_type), intent(out) :: marbl_tracer_metadata(:) ! descriptors for each tracer
- type(marbl_tracer_index_type) , intent(in)  :: marbl_tracer_indices
- type(unit_system_type) , intent(in) :: unit_system
+ type(unit_system_type),    intent(in) :: unit_system
+ type(marbl_tracer_index_type), intent(in) :: marbl_tracer_indices
+ type (marbl_tracer_metadata_type), intent(out) :: marbl_tracer_metadata(:) ! descriptors for each tracer
 
  !-----------------------------------------------------------------------
  ! local variables
@@ -479,15 +479,9 @@ subroutine marbl_init_non_autotroph_tracer_metadata(short_name, long_name, unit_
  marbl_tracer_metadata%long_name = long_name
  if ((trim(short_name) == "ALK") .or. &
  (trim(short_name) == "ALK_ALT_CO2")) then
- if (unit_system%unit_system == 'cgs') then
- marbl_tracer_metadata%units = 'neq/cm^3'
- marbl_tracer_metadata%tend_units = 'neq/cm^3/s'
- marbl_tracer_metadata%flux_units = 'neq/cm^2/s'
- else
- marbl_tracer_metadata%units = 'meq/m^3'
- marbl_tracer_metadata%tend_units = 'meq/m^3/s'
- marbl_tracer_metadata%flux_units = 'meq/m^2/s'
- endif
+ marbl_tracer_metadata%units = unit_system%alk_conc_units
+ marbl_tracer_metadata%tend_units = unit_system%alk_conc_tend_units
+ marbl_tracer_metadata%flux_units = unit_system%alk_conc_flux_units
  else
  marbl_tracer_metadata%units = unit_system%conc_units
  marbl_tracer_metadata%tend_units = unit_system%conc_tend_units
@@ -693,7 +687,7 @@ subroutine marbl_init_surface_flux_forcing_fields(num_elements, surface_flux_for
  if (id .eq. ind%u10_sqr_id) then
  found = .true.
  surface_flux_forcings(id)%metadata%varname = 'u10_sqr'
- write(surface_flux_forcings(id)%metadata%field_units, '(2A)')  trim(unit_system%L), '^2/s^2'
+ write(surface_flux_forcings(id)%metadata%field_units, '(2A)') trim(unit_system%L), '^2/s^2'
  end if
 
  ! Sea-surface salinity

src/marbl_interface.F90

@@ -116,7 +116,7 @@ module marbl_interface
  type(marbl_running_mean_0d_type) , public, allocatable :: glo_scalar_rmean_surface_flux(:)
 
  ! private data
- type(unit_system_type), private  :: unit_system
+ type(unit_system_type), private :: unit_system
  type(marbl_PAR_type), private :: PAR
  type(autotroph_derived_terms_type), private :: autotroph_derived_terms
  type(autotroph_local_type), private :: autotroph_local
@@ -163,6 +163,7 @@ module marbl_interface
  get_string
  procedure, public :: get_settings_var_cnt
  procedure, public :: get_output_for_GCM
+ procedure, public :: get_conc_flux_units
  procedure, private :: inquire_settings_metadata_by_name
  procedure, private :: inquire_settings_metadata_by_id
  procedure, private :: put_real
@@ -199,7 +200,7 @@ subroutine init(this, &
  gcm_delta_z, &
  gcm_zw, &
  gcm_zt, &
- unit_system,  &
+ unit_system_opt, &
  lgcm_has_global_ops)
 
  use marbl_init_mod, only : marbl_init_log_and_timers
@@ -222,11 +223,12 @@ subroutine init(this, &
  real(r8), intent(in) :: gcm_delta_z(gcm_num_levels) ! thickness of layer k
  real(r8), intent(in) :: gcm_zw(gcm_num_levels) ! thickness of layer k
  real(r8), intent(in) :: gcm_zt(gcm_num_levels) ! thickness of layer k
- character(len=*), optional, intent(in) :: unit_system
+ character(len=*), optional, intent(in) :: unit_system_opt
  logical, optional, intent(in) :: lgcm_has_global_ops
 
  character(len=*), parameter :: subname = 'marbl_interface:init'
  integer, parameter :: num_elements_interior_tendency = 1 ! FIXME #66: get this value from interface, let it vary
+ character(len=char_len) :: unit_system_opt_loc
 
  !--------------------------------------------------------------------
  ! initialize status log and timers
@@ -259,7 +261,12 @@ subroutine init(this, &
  ! Initialize parameters that do not depend on tracer count or PFT categories
  !---------------------------------------------------------------------------
 
- call this%unit_system%set(this%StatusLog, unit_system)
+ if (present(unit_system_opt)) then
+ unit_system_opt_loc = unit_system_opt
+ else
+ unit_system_opt_loc = 'cgs'
+ end if
+ call this%unit_system%set(unit_system_opt_loc, this%StatusLog)
  if (this%StatusLog%labort_marbl) then
  call this%StatusLog%log_error_trace("unit_system%set", subname)
  return
@@ -787,7 +794,18 @@ end subroutine get_output_for_GCM
 
  !***********************************************************************
 
- subroutine inquire_settings_metadata_by_name(this, varname, id, lname, units, datatype)
+ function get_conc_flux_units(this)
+
+ class (marbl_interface_class), intent(inout) :: this
+ character(len=char_len) :: get_conc_flux_units
+
+ get_conc_flux_units = this%unit_system%conc_flux_units
+
+ end function get_conc_flux_units
+
+ !***********************************************************************
+
+ subroutine inquire_settings_metadata_by_name(this, varname, id, lname, units, datatype)
 
  class (marbl_interface_class), intent(inout) :: this
  character(len=*), intent(in) :: varname

src/marbl_interface_private_types.F90

@@ -225,20 +225,20 @@ module marbl_interface_private_types
  !***********************************************************************
 
  type, public :: marbl_particulate_share_type
- ! MNL TODO: update units for cgs / mks compatibility
- type(column_sinking_particle_type) :: POC ! base units = nmol C
- type(column_sinking_particle_type) :: POP ! base units = nmol P
- type(column_sinking_particle_type) :: P_CaCO3 ! base units = nmol CaCO3
- type(column_sinking_particle_type) :: P_CaCO3_ALT_CO2 ! base units = nmol CaCO3
- type(column_sinking_particle_type) :: P_SiO2 ! base units = nmol SiO2
- type(column_sinking_particle_type) :: dust ! base units = g
- type(column_sinking_particle_type) :: P_iron ! base units = nmol Fe
+ ! units differ depending for cgs and mks
+ type(column_sinking_particle_type) :: POC ! base units = nmol or mmol C
+ type(column_sinking_particle_type) :: POP ! base units = nmol or mmol P
+ type(column_sinking_particle_type) :: P_CaCO3 ! base units = nmol or mmol CaCO3
+ type(column_sinking_particle_type) :: P_CaCO3_ALT_CO2 ! base units = nmol or mmol CaCO3
+ type(column_sinking_particle_type) :: P_SiO2 ! base units = nmol or mmol SiO2
+ type(column_sinking_particle_type) :: dust ! base units = g or kg
+ type(column_sinking_particle_type) :: P_iron ! base units = nmol or mmol Fe
 
  real(r8), allocatable :: decay_CaCO3_fields (:) ! scaling factor for dissolution of CaCO3
  real(r8), allocatable :: decay_POC_E_fields (:) ! scaling factor for dissolution of excess POC
  real(r8), allocatable :: decay_Hard_fields (:) ! scaling factor for dissolution of Hard Ballast
- real(r8), allocatable :: poc_diss_fields (:) ! diss. length used (cm)
- real(r8), allocatable :: caco3_diss_fields (:) ! caco3 diss. length used (cm)
+ real(r8), allocatable :: poc_diss_fields (:) ! diss. length used (L)
+ real(r8), allocatable :: caco3_diss_fields (:) ! caco3 diss. length used (L)
  real(r8), allocatable :: POC_remin_fields (:) ! POC remin from ecosys before it gets modified for k=KMT
  real(r8), allocatable :: POC_prod_avail_fields (:) ! POC production available for excess POC flux