maecs: improved description of surface O2 flux + fixed dependency chain

d3dcf577 · Onur.Kerimoglu · cb0ed708 · d3dcf577 · d3dcf577 · d3dcf577
Commit d3dcf577 authored Jan 09, 2018 by Onur.Kerimoglu
--- a/src/models/hzg/maecs/maecs.F90
+++ b/src/models/hzg/maecs/maecs.F90
@@ -761,12 +761,13 @@ call self%register_diagnostic_variable(self%id_qualPOM, 'qualPOM','-', 'Quality_
  output=DOUT)
 call self%register_diagnostic_variable(self%id_qualDOM, 'qualDOM','-', 'Quality_of_DOM_ qualDOM', &
  output=DOUT)
-if (self%BioOxyOn) then
-  call self%register_diagnostic_variable(self%id_Denitr,  'Denitr','mmol-N/m**3/d', 'denitrification_rate_ Denitr', output=DOUT)
-  call self%register_diagnostic_variable(self%id_OxCons,  'OxCons','mmol-O2/m**3/d', 'oxygen_consumption_rate_ OxCons', output=DOUT)
-  call self%register_horizontal_diagnostic_variable(self%id_O2flux_diag,'O2flux','mmol-O2/m**2/d','oxygen_flux_between_sea_water_and_air', output=DOUT)
-  call self%register_diagnostic_variable(self%id_no3,     'no3','mmol-N/m**3', 'Nitrate_ no3',output=DOUT)
-end if
+  if (self%BioOxyOn) then
+    call self%register_diagnostic_variable(self%id_Denitr,  'Denitr','mmol-N/m**3/d', 'denitrification_rate_ Denitr', output=DOUT)
+    call self%register_diagnostic_variable(self%id_OxCons,  'OxCons','mmol-O2/m**3/d', 'oxygen_consumption_rate_ OxCons', output=DOUT)
+    call self%register_horizontal_diagnostic_variable(self%id_O2flux_diag,'O2flux','mmol-O2/m**2/d','oxygen_flux_between_sea_water_and_air', output=DOUT)
+    call self%register_diagnostic_variable(self%id_no3,     'no3','mmol-N/m**3', 'Nitrate_ no3',output=DOUT)
+    call self%register_diagnostic_variable(self%id_oxy_percSat,   'O2_percSat','%', 'O2_percent_saturation',output=DOUTa)
+  end if
 end if

 if (self%PhysiolDiagOn) then
@@ -840,6 +841,8 @@ if (self%BGC0DDiagOn .and. .not. self%PhosphorusOn) call self%fatal_error('maecs
 call self%register_dependency(self%id_att,standard_variables%attenuation_coefficient_of_photosynthetic_radiative_flux)
 call self%register_dependency(self%id_attSPM_dep,'attSPM','m-1','SPM_caused_attenuation_in_water')
 call self%register_dependency(self%id_temp,standard_variables%temperature)
+call self%register_dependency(self%id_salt,standard_variables%practical_salinity) !psu
+call self%register_horizontal_dependency(self%id_wind,standard_variables%wind_speed) !m/s
 call self%register_dependency(self%id_par,standard_variables%downwelling_photosynthetic_radiative_flux)
 call self%register_global_dependency(self%id_doy,standard_variables%number_of_days_since_start_of_the_year)
 call self%register_dependency(self%id_vphys_dep,'vphys','','physiological_state_of_sinking_cells')
@@ -863,6 +866,7 @@ if (self%Budget2DDiagOn) then
    call self%register_horizontal_diagnostic_variable(self%id_totS_vertint_diag,'totS_vertint','mmol-Si/m**2','vertical_integral_total_silicate', output=DOUT)
 end if

+!to calculate vertical integrals (somehow doesn't work. check)
 if (self%BGC2DDiagOn) then
    call self%register_dependency(self%id_GPPR_dep,'GPPR','mmol-C/m**3/d','Gross_Primary_Production_Rate')
    call self%register_dependency(self%id_GPPR_vertint,vertical_integral(self%id_GPPR_dep))

--- a/src/models/hzg/maecs/maecs_do.F90
+++ b/src/models/hzg/maecs/maecs_do.F90
@@ -79,7 +79,7 @@ real(rk),parameter :: T0 = 288.15_rk ! reference Temperature fixed to 15 degC
 real(rk),parameter :: Q10b = 1.5_rk
 real(rk) :: Cprod, Nprod, Pprod
 real(rk) :: poc
-real(rk) :: AnoxicMin,Denitrific,OxicMin,Nitri,OduDepo,OduOx,pDepo, Anammox
+real(rk) :: AnoxicMin,Denitrific,OxicMin,Nitri,OduDepo,OduOx,pDepo, Anammox,oxy_sat
 real(rk) :: prodO2, rhochl, uptNH4, uptNO3, uptchl, uptN, respphyto,faeces, min_Cmass
 real(rk) :: vir_max, vir_mu, infect
 logical  :: IsCritical = .false. ! phyC and phyN below reasonable range ?
@@ -142,6 +142,7 @@ end if

 !S_GED
 _GET_(self%id_temp, env%temp)  ! water temperature
+_GET_(self%id_salt, env%salt)  ! salinity
 _GET_(self%id_par, env%par)  ! light photosynthetically active radiation
 if (self%ChemostatOn) then
  if (_AVAILABLE_(self%id_CO2)) then
@@ -748,14 +749,19 @@ if (self%DebugDiagOn) then
 end if
 if (self%BGC0DDiagOn) then
  _SET_DIAGNOSTIC_(self%id_GPPR, _REPLNAN_(phy%gpp*phy%C))   !average gross_primary_production_
-  _SET_DIAGNOSTIC_(self%id_Denitr, _REPLNAN_(0.8*Denitrific)) !average denitrification_rate_
-  _SET_DIAGNOSTIC_(self%id_OxCons, _REPLNAN_(OxCons)) !average denitrification_rate_  
  _SET_DIAGNOSTIC_(self%id_dPAR, _REPLNAN_(env%par))         !average Photosynthetically_Active_Radiation_
  _SET_DIAGNOSTIC_(self%id_DNP, _REPLNAN_(nut%N/(nut%P+self%small))) !average DIN:DIP_ratio_
  _SET_DIAGNOSTIC_(self%id_QNP, _REPLNAN_(phy%Q%N/phy%Q%P))  !average N:P_ratio_
  _SET_DIAGNOSTIC_(self%id_qualPOM, _REPLNAN_(qualPOM))      !average Quality_of_POM_
  _SET_DIAGNOSTIC_(self%id_qualDOM, _REPLNAN_(qualDOM))      !average Quality_of_DOM_
-  _SET_DIAGNOSTIC_(self%id_no3, _REPLNAN_(no3))              !average Nitrate_
+  if (self%BioOxyOn) then
+    _SET_DIAGNOSTIC_(self%id_Denitr, _REPLNAN_(0.8*Denitrific)) !average denitrification_rate_
+    _SET_DIAGNOSTIC_(self%id_OxCons, _REPLNAN_(OxCons)) !average denitrification_rate_
+    _SET_DIAGNOSTIC_(self%id_no3, _REPLNAN_(no3))              !average Nitrate_
+    !calculate %O2 saturation
+    oxy_sat = osat_weiss(env%temp,env%salt)  ! saturation oxygen in mmolO2/m**3 
+    _SET_DIAGNOSTIC_(self%id_oxy_percSat, _REPLNAN_(env%oxy/oxy_sat*100)) !
+  end if
 end if
 if (self%PhysiolDiagOn) then
  _SET_DIAGNOSTIC_(self%id_chl2C, _REPLNAN_(phy%theta*phy%rel_chloropl/12)) !average chlorophyll:carbon_ratio_=_chl-a/chloroplast-C_*_chloroplast-C/phy-molC_*_1molC/12gC_
@@ -936,7 +942,8 @@ subroutine maecs_do_surface(self,_ARGUMENTS_DO_SURFACE_)
   class (type_hzg_maecs), intent(in) :: self
   _DECLARE_ARGUMENTS_DO_SURFACE_

-   real(rk) :: tot_vi_C,tot_vi_N,tot_vi_P,tot_vi_S, O2flux,O2airbl,oxy,tot_vi_GPPR,tot_vi_Denitr,tot_vi_OxCons
+   real(rk) :: tot_vi_C,tot_vi_N,tot_vi_P,tot_vi_S, O2flux,O2airbl,oxy,tot_vi_GPPR,tot_vi_Denitr,tot_vi_OxCons,temp,salt,wnd,sc,p_vel
+   real(rk),parameter  :: secs_pr_day=86400.0_rk
   
 !define _REPLNAN_(X) X !changes back to original code
 #define _REPLNAN_(X) nan_num(X)
@@ -946,7 +953,7 @@ subroutine maecs_do_surface(self,_ARGUMENTS_DO_SURFACE_)
 #if _DEBUG_
 write(*,'(A)') 'begin surface_DO'
 #endif
-
+      
      if (self%BGC2DDiagOn) then
        _GET_HORIZONTAL_(self%id_GPPR_vertint,tot_vi_GPPR)
        _SET_HORIZONTAL_DIAGNOSTIC_(self%id_GPPR_vertint_diag,_REPLNAN_(tot_vi_GPPR))
@@ -981,14 +988,40 @@ write(*,'(A)') 'begin surface_DO'
 ! airsea_ex is the average diffusivity coefficient (m2/sec) divided by the thickness of the boundary layer.
 ! for O2 in mmol m-3, the rate of exchange in mmol m-2 s-1).
 ! Positive values imply a flux into the water, negative: out of the water. 
-         O2airbl = self%O2_sat
-!        _GET_HORIZONTAL_(self%id_O2airbl, O2airbl)! boundary layer dissolved oxygen in mmolO2/m**3
+        
+        !o2 saturation
+        !oxy_sat = self%O2_sat
+        _GET_(self%id_temp, temp)  ! water temperature (C)
+        _GET_(self%id_salt, salt)  ! salinity (psu)
+        oxy_sat = osat_weiss(temp,salt)  ! saturation oxygen in mmolO2/m**3
+        
+        !exchange coefficient (piston velocity) (copied from gotm/ergom)
+        _GET_HORIZONTAL_(self%id_wind,wnd)
+        sc=1450.+(1.1*temp-71.0_rk)*temp
+        if (wnd .gt. 13.0_rk) then
+          p_vel = 5.9_rk*(5.9_rk*wnd-49.3_rk)/sqrt(sc)
+        else
+          if (wnd .lt. 3.6_rk) then
+            p_vel = 1.003_rk*wnd/(sc)**(0.66_rk)
+          else
+            p_vel = 5.9_rk*(2.85_rk*wnd-9.65_rk)/sqrt(sc)
+          end if
+        end if
+        if (p_vel .lt. 0.05_rk) then
+          p_vel = 0.05_rk
+        end if
+        p_vel = p_vel/secs_pr_day
+        
        _GET_(self%id_oxy, oxy)   ! sea water dissolved oxygen in mmolO2/m**3
-
-        O2flux  = self%ex_airsea * (O2airbl - oxy)!
+        !O2flux  = self%ex_airsea * (oxy_sat - oxy)!
+        O2flux =p_vel*(oxy_sat-oxy)
+        !write (*,'(A,5(F10.3))') 'wind,pvel,oxy_sat,O2flux,O2flux_old=', wnd,p_vel,oxy_sat,O2flux,O2flux,self%ex_airsea * (oxy_sat - oxy)
        _SET_SURFACE_EXCHANGE_(self%id_oxy, O2flux )
-        if (self%BGC2DDiagOn) then
+        if (self%BGC0DDiagOn) then
          _SET_HORIZONTAL_DIAGNOSTIC_(self%id_O2flux_diag, _REPLNAN_(O2flux)) ! converts mmol/m2.s to mmol/m2.d
+        end if
+        !calculate vertical integrals
+        if (self%BGC2DDiagOn) then
          _GET_HORIZONTAL_(self%id_Denitr_vertint,tot_vi_Denitr)
          _SET_HORIZONTAL_DIAGNOSTIC_(self%id_Denitr_vertint_diag, _REPLNAN_(tot_vi_Denitr))
          _GET_HORIZONTAL_(self%id_OxCons_vertint,tot_vi_OxCons)

--- a/src/models/hzg/maecs/maecs_functions.F90
+++ b/src/models/hzg/maecs/maecs_functions.F90
@@ -17,7 +17,7 @@
   public   uptflex       ,& 
            queuefunc, queuefunc0, queuederiv ,&
            smooth_small, sinking, min_mass, &
-            calc_sensitivities, calc_internal_states, nan_num
+            calc_sensitivities, calc_internal_states, nan_num, osat_weiss
   
 contains  

@@ -629,7 +629,73 @@ real(rk) function nan_num(x)
   !write(0,*) xnan,x,nan_num
   
   end function nan_num 
+
+!copied from gotm/ergom
+! Weiss formula for the saturation oxygen (osat)
+   real(rk) function osat_weiss(t,s)
+!
+! !DESCRIPTION:
+! Weiss formula for the saturation oxygen (osat) \cite{Weiss1970}:
+!
+! \begin{equation}\label{osat_weiss}
+! O_{sat}= \exp\left[a_1 +a_2\frac{100}{T}+a_3\ln\left(\frac{T}{100}\right)
+! +a_4\frac{T}{100}+S \left\{b_1+b_2\frac{T}{100}
+! +b_3\left(\frac{T}{100}\right)^2   \right\}\right],
+! \end{equation}
+!
+! where $T$ is the temperature in Kelvin and the empirical constants are
+! given in table \ref{table_weiss}.
+!
+! \begin{table}[h]
+! \begin{center}
+! \begin{tabular}{|l|l|l|l|l|l|l|}
+! \hline
+! $a_1$ & $a_2$ & $a_3$ & $a_4$ & $b_1$ & $b_2$ & $b_3$ \\ \hline
+! -173.4292 &
+! 249.6339 &
+! 143.3483 &
+! -21.8492 &
+! -0.033096 &
+! 0.014259 &
+! -0.001700 \\ \hline
+! \end{tabular}
+! \caption{Constants for the oxygen saturation formula by \cite{Weiss1970},
+! see equation (\ref{osat_weiss}).}
+! \label{table_weiss}
+! \end{center}
+! \end{table}
+!
+! !USES:
+   IMPLICIT NONE
+!
+! !INPUT PARAMETERS:
+!  type (type_gotm_ergom), intent(in) :: self
+  real(rk), intent(in)                 :: t,s
+!
+! !REVISION HISTORY:
+!  Original author(s): Hans Burchard, Karsten Bolding
+!
+! !LOCAL VARIABLES:
+  real(rk)                 :: tk
+  real(rk)                 :: aa1=-173.4292_rk
+  real(rk)                 :: aa2=249.6339_rk
+  real(rk)                 :: a3=143.3483_rk
+  real(rk)                 :: a4=-21.8492_rk
+  real(rk)                 :: b1=-0.033096_rk
+  real(rk)                 :: b2=0.014259_rk
+  real(rk)                 :: b3=-0.001700_rk
+  real(rk)                 :: kelvin=273.16_rk
+  real(rk)                 :: mol_per_liter=44.661_rk
+!EOP
+!-----------------------------------------------------------------------
+!BOC
+   tk=(t+kelvin)*0.01_rk
+   osat_weiss=exp(aa1+aa2/tk+a3*log(tk)+a4*tk    &
+              +s*(b1+(b2+b3*tk)*tk))*mol_per_liter
+   return
+   end function osat_weiss
   
+!EOC
 end module maecs_functions
 !------------------------------------------------------

--- a/src/models/hzg/maecs/maecs_types.F90
+++ b/src/models/hzg/maecs/maecs_types.F90
@@ -36,10 +36,12 @@ type (type_maecs_nutindex) :: nutind
 ! standard fabm model types
 type (type_state_variable_id)        :: id_nutN,id_nutP,id_nutS,id_phyC,id_phyN,id_phyP,id_phyS,id_zooC,id_detC,id_detN,id_detP,id_detS,id_domC,id_domN,id_domP,id_RNit,id_Rub,id_chl,id_nh3,id_oxy,id_odu
 type (type_dependency_id)            :: id_temp
+type (type_dependency_id)            :: id_salt
 type (type_dependency_id)            :: id_par
 type (type_dependency_id)            :: id_att
 type (type_dependency_id)            :: id_CO2
-type (type_global_dependency_id)            :: id_doy
+type (type_global_dependency_id)     :: id_doy
+type (type_horizontal_dependency_id) :: id_wind
 type (type_dependency_id)            :: id_totC
 type (type_horizontal_dependency_id)            :: id_totC_vertint
 type (type_dependency_id)            :: id_totN
@@ -62,7 +64,7 @@ type (type_horizontal_diagnostic_variable_id)            :: id_totP_vertint_diag
 type (type_horizontal_diagnostic_variable_id)            :: id_totS_vertint_diag
 type (type_horizontal_dependency_id)            :: id_zmax
 type (type_horizontal_diagnostic_variable_id)            :: id_O2flux_diag
-type (type_diagnostic_variable_id)   :: id_datt, id_dattSPM, id_vphys, id_GPPR, id_Denitr, id_OxCons, id_dPAR, id_chl2C, id_Theta, id_fracR, id_fracT, id_fracNU, id_DNP, id_QNP, id_QN, id_QP, id_QSi, id_aVN, id_aVP, id_aVSi, id_faN, id_faP, id_faSi, id_rQN, id_rQP, id_rQSi, id_tmp, id_fac1, id_fac2, id_fac3, id_fac4, id_fac5, id_phyUR, id_phyRER, id_phyELR, id_phyALR, id_phyVLR, id_phyGLR, id_vsinkr, id_zoomort, id_qualPOM, id_qualDOM, id_no3, id_UCpot
+type (type_diagnostic_variable_id)   :: id_datt, id_dattSPM, id_vphys, id_GPPR, id_Denitr, id_OxCons, id_dPAR, id_chl2C, id_Theta, id_fracR, id_fracT, id_fracNU, id_DNP, id_QNP, id_QN, id_QP, id_QSi, id_aVN, id_aVP, id_aVSi, id_faN, id_faP, id_faSi, id_rQN, id_rQP, id_rQSi, id_tmp, id_fac1, id_fac2, id_fac3, id_fac4, id_fac5, id_phyUR, id_phyRER, id_phyELR, id_phyALR, id_phyVLR, id_phyGLR, id_vsinkr, id_zoomort, id_qualPOM, id_qualDOM, id_no3, id_UCpot, id_oxy_percSat
 real(rk) ::  nutN_initial, nutP_initial, nutS_initial, phyC_initial, phyN_initial, phyP_initial, phyS_initial, zooC_initial, detC_initial, detN_initial, detP_initial, detS_initial, domC_initial, domN_initial, domP_initial, RNit_initial, frac_Rub_ini, frac_chl_ini, nh3_initial, oxy_initial, odu_initial
 real(rk) ::  P_max, alpha, sigma, theta_LHC, rel_chloropl_min, QN_phy_0, QN_phy_max, V_NC_max, AffN, zeta_CN, zstoich_PN, exud_phy, QP_phy_0, QP_phy_max, V_PC_max, AffP, QSi_phy_0, QSi_phy_max, V_SiC_max, AffSi, MaxRelQ, syn_nut, adap_rub, adap_theta, tau_regV, disease, mort_ODU, decay_pigm, decay_nut, phi_agg, agg_doc, vir_loss, vir_bmass, sink_phys, vS_phy, vS_phy0, vS_det, hydrol, remin, Nqual, remNP, denit, PON_denit, Q10, T_ref, NutOrder
 real(rk) ::  const_NC_zoo, const_PC_zoo, g_max, k_grazC, yield_zoo, basal_resp_zoo, mort_zoo, zm_fa_delmax, zm_fa_inf, Q10z, fT_exp_mort
@@ -75,7 +77,7 @@ end type type_maecs_base_model

 !
 type type_maecs_env
- real(rk) :: RNit, nh3, oxy, odu,  temp,par,CO2,doy,attf_dep,vphys_dep,GPPR_dep,GPPR_vertint,GPPR_vertint_diag,Denitr_dep,Denitr_vertint,Denitr_vertint_diag,OxCons_dep,OxCons_vertint,Oxcons_vertint_diag,zmax,O2flux_diag
+ real(rk) :: RNit, nh3, oxy, odu,  temp,salt,par,CO2,doy,attf_dep,vphys_dep,GPPR_dep,GPPR_vertint,GPPR_vertint_diag,Denitr_dep,Denitr_vertint,Denitr_vertint_diag,OxCons_dep,OxCons_vertint,Oxcons_vertint_diag,zmax,O2flux_diag
 end type
 type type_maecs_rhs
 real(rk) :: nutN,nutP,nutS,phyC,phyN,phyP,phyS,zooC,detC,detN,detP,detS,domC,domN,domP,RNit,Rub,chl,nh3,oxy,odu