Integrate sed* tracers over configurable bottom thickness (#42)

Co-authored-by: Pearse J Buchanan <[email protected]>

Author: dougiesquire

generic_tracers/generic_WOMBATlite.F90

@@ -165,6 +165,7 @@ module generic_WOMBATlite
         zoolmor, &
         zooqmor, &
         detlrem, &
+        bottom_thickness, &
         detlrem_sed, &
         wdetbio, &
         wdetmax, &
@@ -1148,7 +1149,7 @@ subroutine generic_WOMBATlite_register_diag(diag_list)
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'seddep', 'Depth of the sediment', 'h', '1', 's', 'm', 'f')
+        'seddep', 'Depth of the bottom layer', 'h', '1', 's', 'm', 'f')
     wombat%id_seddep = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
@@ -1158,42 +1159,42 @@ subroutine generic_WOMBATlite_register_diag(diag_list)
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'sedtemp', 'Temperature at the deepest grid cell', 'h', '1', 's', 'deg C', 'f')
+        'sedtemp', 'Temperature in the bottom layer', 'h', '1', 's', 'deg C', 'f')
     wombat%id_sedtemp = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'sedsalt', 'Salinity at the deepest grid cell', 'h', '1', 's', 'psu', 'f')
+        'sedsalt', 'Salinity in the bottom layer', 'h', '1', 's', 'psu', 'f')
     wombat%id_sedsalt = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'sedno3', 'Nitrate at the deepest grid cell', 'h', '1', 's', 'mol/kg', 'f')
+        'sedno3', 'Nitrate concentration in the bottom layer', 'h', '1', 's', 'mol/kg', 'f')
     wombat%id_sedno3 = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'seddic', 'Dissolved inorganic carbon at the deepest grid cell', 'h', '1', 's', 'mol/kg', 'f')
+        'seddic', 'Dissolved inorganic carbon concentration in the bottom layer', 'h', '1', 's', 'mol/kg', 'f')
     wombat%id_seddic = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'sedalk', 'Alkalinity at the deepest grid cell', 'h', '1', 's', 'mol/kg', 'f')
+        'sedalk', 'Alkalinity concentration in the bottom layer', 'h', '1', 's', 'mol/kg', 'f')
     wombat%id_sedalk = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'sedhtotal', 'H+ ions at the deepest grid cell', 'h', '1', 's', 'mol/kg', 'f')
+        'sedhtotal', 'H+ ion concentration in the bottom layer', 'h', '1', 's', 'mol/kg', 'f')
     wombat%id_sedhtotal = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'sedco3', 'CO3 ions at the deepest grid cell', 'h', '1', 's', 'mol/kg', 'f')
+        'sedco3', 'CO3 ion concentration in the bottom layer', 'h', '1', 's', 'mol/kg', 'f')
     wombat%id_sedco3 = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
     vardesc_temp = vardesc( &
-        'sedomega_cal', 'Calcite saturation state at the deepest grid cell', 'h', '1', 's', ' ', 'f')
+        'sedomega_cal', 'Calcite saturation state in the bottom layer', 'h', '1', 's', ' ', 'f')
     wombat%id_sedomega_cal = register_diag_field(package_name, vardesc_temp%name, axes(1:2), &
         init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
 
@@ -1418,7 +1419,12 @@ subroutine user_add_params
     ! Detritus maximum sinking rate coefficient [m/s]
     !-----------------------------------------------------------------------
     call g_tracer_add_param('wdetmax', wombat%wdetmax, 36.0/86400.0)
-    
+
+    ! Bottom thickness [m]
+    !-----------------------------------------------------------------------
+    ! Thickness over which tracer values are integrated to define the bottom layer
+    call g_tracer_add_param('bottom_thickness', wombat%bottom_thickness, 1.0)
+
     ! Detritus remineralisation rate constant in sediments [1/s]
     !-----------------------------------------------------------------------
     ! This would normally equal detlrem, but we set the default value to be
@@ -2010,7 +2016,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt,  &
     real, dimension(:,ilb:,jlb:,:), intent(in) :: opacity_band
 
     integer                                 :: isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau, tn
-    integer                                 :: i, j, k, n, nz
+    integer                                 :: i, j, k, n, nz, k_bot
     real, dimension(:,:,:), pointer         :: grid_tmask
     integer, dimension(:,:), pointer        :: grid_kmt
     integer, dimension(:,:), allocatable    :: kmeuph ! deepest level of euphotic zone
@@ -2044,6 +2050,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt,  &
     real                                    :: zoo_slmor, epsmin
     real                                    :: hco3, diss_cal, diss_ara, diss_det
     real                                    :: avedetbury, avecaco3bury
+    real                                    :: dzt_bot, dzt_bot_os
     real, dimension(:,:,:,:), allocatable   :: n_pools, c_pools
     logical                                 :: used
 
@@ -3006,6 +3013,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt,  &
           print *, "       Nested timestep number =", tn
           print *, " "
           print *, "       Biological N budget (molN/kg) at two timesteps =", n_pools(i,j,k,1), n_pools(i,j,k,2)
+          print *, "       Difference in budget between timesteps =", n_pools(i,j,k,2) - n_pools(i,j,k,1)
           print *, " "
           print *, "       NO3 (molNO3/kg) =", wombat%f_no3(i,j,k)
           print *, "       PHY (molN/kg) =", wombat%f_phy(i,j,k) * 16.0 / 122.0
@@ -3030,6 +3038,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt,  &
           print *, "       Nested timestep number =", tn
           print *, " "
           print *, "       Biological C budget (molC/kg) at two timesteps =", c_pools(i,j,k,1), c_pools(i,j,k,2)
+          print *, "       Difference in budget between timesteps =", c_pools(i,j,k,2) - c_pools(i,j,k,1)
           print *, " "
           print *, "       DIC (molC/kg) =", wombat%f_dic(i,j,k)
           print *, "       ALK (molC/kg) =", wombat%f_alk(i,j,k)
@@ -3188,26 +3197,65 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt,  &
     call g_tracer_get_pointer(tracer_list, 'detfe_sediment', 'field', wombat%p_detfe_sediment) ! [mol/m2]
     call g_tracer_get_pointer(tracer_list, 'caco3_sediment', 'field', wombat%p_caco3_sediment) ! [mol/m2]
 
+    ! Get bottom conditions, including those that influence bottom fluxes. Bottom conditions are
+    ! calculated over a layer defined by wombat%bottom_thickness (default 1 m). This is done because
+    ! the bottom layers in MOM6 are usually "vanished" layers. This approach is based on what is done
+    ! in COBALT v3.
     do j = jsc,jec; do i = isc,iec;
+      k_bot = 0
+      dzt_bot = 0.0
+      do k = grid_kmt(i,j),1,-1
+        if (dzt_bot .lt. wombat%bottom_thickness) then
+          k_bot = k
+          dzt_bot = dzt_bot + dzt(i,j,k) ! [m]
+          wombat%sedtemp(i,j) = wombat%sedtemp(i,j) + Temp(i,j,k) * dzt(i,j,k) ! [m*degC]
+          wombat%sedsalt(i,j) = wombat%sedsalt(i,j) + Salt(i,j,k) * dzt(i,j,k) ! [m*psu]
+          wombat%sedno3(i,j) = wombat%sedno3(i,j) + wombat%f_no3(i,j,k) * dzt(i,j,k) ! [m*mol/kg]
+          wombat%seddic(i,j) = wombat%seddic(i,j) + wombat%f_dic(i,j,k) * dzt(i,j,k) ! [m*mol/kg]
+          wombat%sedalk(i,j) = wombat%sedalk(i,j) + wombat%f_alk(i,j,k) * dzt(i,j,k) ! [m*mol/kg]
+          wombat%sedhtotal(i,j) = wombat%sedhtotal(i,j) + wombat%htotal(i,j,k) * dzt(i,j,k) ! [m*mol/kg]
+        endif
+      enddo
+      ! Subtract off overshoot
+      dzt_bot_os = dzt_bot - wombat%bottom_thickness
+      wombat%sedtemp(i,j) = wombat%sedtemp(i,j) - Temp(i,j,k_bot) * dzt_bot_os ! [m*degC]
+      wombat%sedsalt(i,j) = wombat%sedsalt(i,j) - Salt(i,j,k_bot) * dzt_bot_os ! [m*psu]
+      wombat%sedno3(i,j) = wombat%sedno3(i,j) - wombat%f_no3(i,j,k_bot) * dzt_bot_os ! [m*mol/kg]
+      wombat%seddic(i,j) = wombat%seddic(i,j) - wombat%f_dic(i,j,k_bot) * dzt_bot_os ! [m*mol/kg]
+      wombat%sedalk(i,j) = wombat%sedalk(i,j) - wombat%f_alk(i,j,k_bot) * dzt_bot_os ! [m*mol/kg]
+      wombat%sedhtotal(i,j) = wombat%sedhtotal(i,j) - wombat%htotal(i,j,k_bot) * dzt_bot_os ! [m*mol/kg]
+      ! Convert to mol/kg
+      wombat%sedtemp(i,j) = wombat%sedtemp(i,j) / wombat%bottom_thickness ! [degC]
+      wombat%sedsalt(i,j) = wombat%sedsalt(i,j) / wombat%bottom_thickness ! [psu]
+      wombat%sedno3(i,j) = wombat%sedno3(i,j) / wombat%bottom_thickness ! [mol/kg]
+      wombat%seddic(i,j) = wombat%seddic(i,j) / wombat%bottom_thickness ! [mol/kg]
+      wombat%sedalk(i,j) = wombat%sedalk(i,j) / wombat%bottom_thickness ! [mol/kg]
+      wombat%sedhtotal(i,j) = wombat%sedhtotal(i,j) / wombat%bottom_thickness ! [mol/kg]
+
+      ! Set seddep as full depth minus half the bottom thickness and sedmask from bottom layer
       k = grid_kmt(i,j)
       if (k .gt. 0) then
-      ! Get bottom conditions: mask, PO4, temp, salt, DIC, Alk, H+ ion
-        wombat%seddep(i,j) = wombat%zw(i,j,k)
+        wombat%seddep(i,j) = max(0.0, wombat%zw(i,j,k) - (wombat%bottom_thickness / 2.0))
         wombat%sedmask(i,j) = grid_tmask(i,j,k)
-        wombat%sedtemp(i,j) = Temp(i,j,k)
-        wombat%sedsalt(i,j) = Salt(i,j,k)
-        wombat%sedno3(i,j) = wombat%f_no3(i,j,k)
-        wombat%seddic(i,j) = wombat%f_dic(i,j,k) + wombat%p_det_sediment(i,j,1) / wombat%Rho_0  ![mol/kg] 
-        wombat%sedalk(i,j) = wombat%f_alk(i,j,k)
-        wombat%sedhtotal(i,j) = wombat%htotal(i,j,k)
       endif
+
+      ! pjb: Sum the water column concentration of DIC and the organic carbon content of the
+      ! sediment to approximate the interstitial (i.e., porewater) DIC concentration.
+      ! We assume that the organic carbon content of the sediment (p_det_sediment) in mol/m2 is
+      ! relevant over one meter, and therefore can be automatically converted to mol/m3 and then
+      ! subsequently converted through the mol/kg using Rho_0. With this assumption these arrays
+      ! can be added together.
+      ! We add these arrays together to simulate the reducing conditions of organic-rich sediments,
+      ! and to calculate a lower omega for calcite, which ensures greater rates of dissolution of
+      ! CaCO3 within the sediment as organic matter accumulates.
+      wombat%seddic(i,j) = wombat%seddic(i,j) + wombat%p_det_sediment(i,j,1) / wombat%Rho_0
     enddo; enddo
 
     call FMS_ocmip2_co2calc(CO2_dope_vec, wombat%sedmask(:,:), &
         wombat%sedtemp(:,:), wombat%sedsalt(:,:), &
         min(wombat%dic_max*mmol_m3_to_mol_kg, max(wombat%seddic(:,:), wombat%dic_min*mmol_m3_to_mol_kg)), &
         max(wombat%sedno3(:,:) / 16., 1e-9), &
-        wombat%sio2(:,:), &
+        wombat%sio2(:,:), & ! dts: This is currently constant, equal to wombat%sio2_surf
         min(wombat%alk_max*mmol_m3_to_mol_kg, max(wombat%sedalk(:,:), wombat%alk_min*mmol_m3_to_mol_kg)), &
         wombat%sedhtotal(:,:)*wombat%htotal_scale_lo, &
         wombat%sedhtotal(:,:)*wombat%htotal_scale_hi, &