diff --git a/source/model/init_var_mod.f90 b/source/model/init_var_mod.f90
index c0d838cd48e12d4a7df0385c161d08b11cc24bb9..09f7083a608dd1bcbc69019e63787453607efc35 100644
--- a/source/model/init_var_mod.f90
+++ b/source/model/init_var_mod.f90
@@ -289,10 +289,10 @@
& Henry_temp_ref, Tb0 + Kelvin0) ! Saturated concentration at the bottom
do i = 1, M+1
if (i <= i_ML) then
- qwater(i) = ch4_atm0 ! Atmospheric concentration
+ qwater(i) = 0. !ch4_atm0 ! Atmospheric concentration
else
- qwater(i) = qwater(max(i_ML,1)) + &
- & (ch4_bot - qwater(max(i_ML,1)))*(h_ML0 - dzeta_int(i)*h10)/(h_ML0 - h10)
+ qwater(i) = 0. !qwater(max(i_ML,1)) + &
+ !& (ch4_bot - qwater(max(i_ML,1)))*(h_ML0 - dzeta_int(i)*h10)/(h_ML0 - h10)
endif
enddo
else
@@ -351,9 +351,9 @@
if (init_T /= 3) then
do i = 1, M+1
if (i <= i_ML) then
- oxyg(i) = o2_atm0 ! Atmospheric concentration
+ oxyg(i) = 0. !o2_atm0 ! Atmospheric concentration
else
- oxyg(i) = oxyg(max(i_ML,1)) - oxyg(max(i_ML,1))*(h_ML0 - dzeta_int(i)*h10)/(h_ML0 - h10)
+ oxyg(i) = 0. !oxyg(max(i_ML,1)) - oxyg(max(i_ML,1))*(h_ML0 - dzeta_int(i)*h10)/(h_ML0 - h10)
endif
oxyg(i) = max(oxyg(i),small_value)
enddo
diff --git a/source/model/lake.f90 b/source/model/lake.f90
index 9bce7427feba6c3b5e47c27b5f548f0050d2b2c5..b4c638dcaaad7b85a61b08480d5bc781f6adecab 100644
--- a/source/model/lake.f90
+++ b/source/model/lake.f90
@@ -989,20 +989,6 @@ if (soilswitch%par == 1) then
& bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso)
endif
- gs%isoilcol = nsoilcols
- call METHANE &
- & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
- & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
- & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
- & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
- & lsm, bathymwater, &
- & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
- & plant_sum,bull_sum,oxid_sum,rprod_sum, &
- & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
- & ice_meth_oxid_total, &
- & h_talik,tot_ice_meth_bubbles,add_to_winter)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
!call METHANE2 & ! two-meth
!& (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
!& ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
@@ -1012,11 +998,22 @@ if (soilswitch%par == 1) then
!& h_talik,tot_ice_meth_bubbles2) ! two-meth
!qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
!qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
-else
- qwater(:,2) = qwater(:,1)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
endif
+gs%isoilcol = nsoilcols
+call METHANE &
+& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
+& fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
+& rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
+& ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
+& lsm, bathymwater, &
+& fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
+& plant_sum,bull_sum,oxid_sum,rprod_sum, &
+& anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
+& ice_meth_oxid_total, &
+& h_talik,tot_ice_meth_bubbles,add_to_winter)
+qmethane(1:M+1) = qwater(1:M+1,2)
+qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
+
!call ADDOXPRODCONS(M,i_maxN,H_mixed_layer,H_photic_zone,prodox,resp,bod,sod,RadWater,dt,ls,gas)
call OXYGEN (gs, Tw2, fbbleflx_o2_sum, fbbleflx_o2(0,nsoilcols), lamoxyg, gsp, fracb0, pressure, wind10, &
& o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2) !Foxyg1) !Foxyg1 added by DG
@@ -1368,20 +1365,6 @@ if (flag_snow == 1) then
& bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso)
endif
- gs%isoilcol = nsoilcols
- call METHANE &
- & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
- & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
- & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
- & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false. ,&
- & lsm, bathymwater, &
- & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
- & plant_sum,bull_sum,oxid_sum,rprod_sum, &
- & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
- & ice_meth_oxid_total, &
- & h_talik,tot_ice_meth_bubbles,add_to_winter)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
! call METHANE2 & ! two-meth
! & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
! & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
@@ -1391,11 +1374,21 @@ if (flag_snow == 1) then
! & h_talik,tot_ice_meth_bubbles2) ! two-meth
! qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
! qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
- else
- qwater(:,2) = qwater(:,1)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
endif
+ gs%isoilcol = nsoilcols
+ call METHANE &
+ & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
+ & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
+ & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
+ & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false. ,&
+ & lsm, bathymwater, &
+ & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
+ & plant_sum,bull_sum,oxid_sum,rprod_sum, &
+ & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
+ & ice_meth_oxid_total, &
+ & h_talik,tot_ice_meth_bubbles,add_to_winter)
+ qmethane(1:M+1) = qwater(1:M+1,2)
+ qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
!call ADDOXPRODCONS(M,i_maxN,H_mixed_layer,H_photic_zone,prodox,resp,bod,sod,RadWater,dt,ls,gas)
call OXYGEN (gs, Tw2, fbbleflx_o2_sum, fbbleflx_o2(0,nsoilcols), lamoxyg, gsp, fracb0, pressure, wind10, &
& o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2) !, Foxyg1) !Foxyg1 added by DG
@@ -1632,20 +1625,6 @@ else
& bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso)
endif
- gs%isoilcol = nsoilcols
- call METHANE &
- & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
- & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
- & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
- & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
- & lsm, bathymwater, &
- & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
- & plant_sum,bull_sum,oxid_sum,rprod_sum, &
- & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
- & ice_meth_oxid_total, &
- & h_talik,tot_ice_meth_bubbles,add_to_winter)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
! call METHANE2 & ! two-meth
! & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
! & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
@@ -1655,11 +1634,22 @@ else
! & h_talik,tot_ice_meth_bubbles2) ! two-meth
! qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
! qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
- else
- qwater(:,2) = qwater(:,1)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
endif
+ gs%isoilcol = nsoilcols
+ call METHANE &
+ & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
+ & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
+ & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
+ & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
+ & lsm, bathymwater, &
+ & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
+ & plant_sum,bull_sum,oxid_sum,rprod_sum, &
+ & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
+ & ice_meth_oxid_total, &
+ & h_talik,tot_ice_meth_bubbles,add_to_winter)
+ qmethane(1:M+1) = qwater(1:M+1,2)
+ qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
+
!call ADDOXPRODCONS(M,i_maxN,H_mixed_layer,H_photic_zone,prodox,resp,bod,sod,RadWater,dt,ls,gas)
call OXYGEN (gs, Tw2, fbbleflx_o2_sum, fbbleflx_o2(0,nsoilcols), lamoxyg, gsp, fracb0, pressure, wind10, &
& o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2) !, Foxyg1) !Foxyg1 added by DG
@@ -1829,20 +1819,6 @@ if (flag_snow == 1) then
& snowmass_init,a,b,c,d,Temp,phi,fetch,dt)
if (soilswitch%par == 1) then
call BUBBLE_BLOCK
- gs%isoilcol = nsoilcols
- call METHANE &
- & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
- & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
- & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
- & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
- & lsm, bathymwater, &
- & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
- & plant_sum,bull_sum,oxid_sum,rprod_sum, &
- & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
- & ice_meth_oxid_total, &
- & h_talik,tot_ice_meth_bubbles,add_to_winter)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
! call METHANE2 & ! two-meth
! & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
! & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
@@ -1853,6 +1829,21 @@ if (flag_snow == 1) then
! qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
! qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
endif
+ gs%isoilcol = nsoilcols
+ call METHANE &
+ & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
+ & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
+ & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
+ & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
+ & lsm, bathymwater, &
+ & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
+ & plant_sum,bull_sum,oxid_sum,rprod_sum, &
+ & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
+ & ice_meth_oxid_total, &
+ & h_talik,tot_ice_meth_bubbles,add_to_winter)
+ qmethane(1:M+1) = qwater(1:M+1,2)
+ qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
+
l2 = l1 + dhi
h2 = h1 + dhw
@@ -1910,20 +1901,6 @@ else
if (soilswitch%par == 1) then
call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
call BUBBLE_BLOCK
- gs%isoilcol = nsoilcols
- call METHANE &
- & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
- & fbbleflx_ch4(0,nsoilcols), wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
- & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
- & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
- & lsm, bathymwater, &
- & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
- & plant_sum,bull_sum,oxid_sum,rprod_sum, &
- & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
- & ice_meth_oxid_total, &
- & h_talik,tot_ice_meth_bubbles,add_to_winter)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
! call METHANE2 & ! two-meth
! & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
! & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
@@ -1934,7 +1911,21 @@ else
! qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
! qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
endif
-
+ gs%isoilcol = nsoilcols
+ call METHANE &
+ & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
+ & fbbleflx_ch4(0,nsoilcols), wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
+ & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
+ & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
+ & lsm, bathymwater, &
+ & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
+ & plant_sum,bull_sum,oxid_sum,rprod_sum, &
+ & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
+ & ice_meth_oxid_total, &
+ & h_talik,tot_ice_meth_bubbles,add_to_winter)
+ qmethane(1:M+1) = qwater(1:M+1,2)
+ qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
+
l2 = l1 + dhi
h2 = h1 + dhw
ls2 = 0.
@@ -1969,20 +1960,6 @@ if4: IF (layer_case==4) THEN
& snowmass_init,a,b,c,d,Temp,phi,fetch,dt)
if (soilswitch%par == 1) then
call BUBBLE_BLOCK
- gs%isoilcol = nsoilcols
- call METHANE &
- & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
- & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
- & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
- & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
- & lsm, bathymwater, &
- & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
- & plant_sum,bull_sum,oxid_sum,rprod_sum, &
- & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
- & ice_meth_oxid_total, &
- & h_talik,tot_ice_meth_bubbles,add_to_winter)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
! call METHANE2 & ! two-meth
! & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
! & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
@@ -1993,6 +1970,21 @@ if4: IF (layer_case==4) THEN
! qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
! qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
endif
+ gs%isoilcol = nsoilcols
+ call METHANE &
+ & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
+ & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
+ & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
+ & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
+ & lsm, bathymwater, &
+ & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
+ & plant_sum,bull_sum,oxid_sum,rprod_sum, &
+ & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
+ & ice_meth_oxid_total, &
+ & h_talik,tot_ice_meth_bubbles,add_to_winter)
+ qmethane(1:M+1) = qwater(1:M+1,2)
+ qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
+
else
ice = 0; snow = 0; water = 0; deepice = 0
@@ -2008,19 +2000,6 @@ if4: IF (layer_case==4) THEN
if (soilswitch%par == 1) then
call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
call BUBBLE_BLOCK
- gs%isoilcol = nsoilcols
- call METHANE(gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
- & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
- & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
- & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
- & lsm, bathymwater, &
- & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
- & plant_sum,bull_sum,oxid_sum,rprod_sum, &
- & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
- & ice_meth_oxid_total, &
- & h_talik,tot_ice_meth_bubbles,add_to_winter)
- qmethane(1:M+1) = qwater(1:M+1,2)
- qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
! call METHANE2 & ! two-meth
! & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
! & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
@@ -2031,6 +2010,20 @@ if4: IF (layer_case==4) THEN
! qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
! qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
endif
+ gs%isoilcol = nsoilcols
+ call METHANE(gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
+ & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
+ & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
+ & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
+ & lsm, bathymwater, &
+ & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
+ & plant_sum,bull_sum,oxid_sum,rprod_sum, &
+ & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
+ & ice_meth_oxid_total, &
+ & h_talik,tot_ice_meth_bubbles,add_to_winter)
+ qmethane(1:M+1) = qwater(1:M+1,2)
+ qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
+
endif
ENDIF if4
diff --git a/source/model/methane_mod.f90 b/source/model/methane_mod.f90
index e3a9c6355e1d2dbcecabbec537c629257804b0c4..8861b965ca71e10f991ca75ca6472eee70d27c45 100644
--- a/source/model/methane_mod.f90
+++ b/source/model/methane_mod.f90
@@ -22,7 +22,8 @@ use DRIVING_PARAMS , only : &
& missing_value, &
& thermokarst_meth_prod, & ! Switch for old organics methane production for thermokarst lakes
& tricemethhydr, &
-& soil_meth_prod ! Switch for methane production from new organics
+& soil_meth_prod, & ! Switch for methane production from new organics
+& soilswitch
use PHYS_CONSTANTS, only : &
& row0, g, row0, roa0, roi, &
@@ -227,6 +228,12 @@ implicit none
real(kind=ireals), pointer :: qsoil(:), qwater(:,:), Twater(:), u(:), v(:)
+
+ if (soilswitch%par == 0 .and. botflux_ch4 == missing_value) then
+ print*, 'If sediments are deactivated, methane flux at bottom must be prescribed: STOP'
+ stop
+ endif
+
! Constant of methane production in sediments is set differently for
! sediments in photic zone and underneath
dz_ = ls%H_photic_zone - h1
@@ -248,6 +255,8 @@ implicit none
! dhw = 0.; dhw0 = 0.
+ soilswitch_if: if (soilswitch%par == 1) then
+
z(:) = zsoil(:)
roots(:) = rootss(:)
@@ -601,6 +610,55 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.)
q_old(i) = qsoil(i)
end do
+! ============================================================
+! the methane flux due to plant-mediated transport
+! =================fplant(t)==================================
+! do i = ms,ml
+! qplant(i) = rkp*tveg*froot(i)*fgrow*
+! * (1.-0.5)*qsoil(i)
+! end do
+
+!=============================================================
+! the methane flux due to ebullition
+!=================febul(t)====================================
+
+ qebul(:) = 0.
+ do i = 1, gs%ns !ms, ml
+ if (wl(i) .gt. anox_crit(i)) then
+ qebul(i) = bull(i)*(q_old(i) - ch4_crit(i)) ! cthresh -> ch4_crit(i), according to explicit scheme
+ endif
+ end do
+
+!Methane bubble flux at the bottom
+ febul0 = 0.
+ do i = 2, gs%ns-1 !ms+1, ml
+! if(wl(i) .gt. anox_crit(i)) then
+! if(z(i).ge.w) then
+! The original statement is replaced by the grid-consistent one
+! febul = febul + 0.5*(z(i)-z(i-1))*(qebul(i)+qebul(i-1))
+ febul0 = febul0 + 0.5*(z(i+1)-z(i-1))*qebul(i)
+! end if
+ end do
+ febul0 = febul0 + qebul(1)*0.5*(z(2)-z(1))
+ febul0 = febul0 + qebul(gs%ns)*0.5*(z(gs%ns)-z(gs%ns-1))
+
+ else
+
+ febul0 = 0.
+
+ endif soilswitch_if
+
+ if (deepestsoil) then
+ ! Adding bubble flux from the lowest soil column to horizontally averaged bubble flux
+ if (gs%nsoilcols > 1) then
+ call BUBBLEFLUXAVER(gs%ix,gs%iy,gs%nsoilcols)
+ else
+ fbbleflx_ch4_sum(0) = febul0*fbbleflx_ch4(0) *bathymwater(1) %area_int
+ fbbleflx_ch4_sum(1:gs%M) = febul0*fbbleflx_ch4(1:gs%M)*bathymwater(1:gs%M)%area_half
+ fbbleflx_ch4_sum(gs%M+1) = febul0*fbbleflx_ch4(gs%M+1)*bathymwater(gs%M+1)%area_int
+ endif
+ endif
+
! =============================================================
! set up the tridiagonal system of equations
! =============================================================
@@ -885,53 +943,24 @@ endif
!print*, Flux_atm
-! ============================================================
-! the methane flux due to plant-mediated transport
-! =================fplant(t)==================================
-! do i = ms,ml
-! qplant(i) = rkp*tveg*froot(i)*fgrow*
-! * (1.-0.5)*qsoil(i)
-! end do
-
-!=============================================================
-! the methane flux due to ebullition
-!=================febul(t)====================================
-
- qebul(:) = 0.
- do i = 1, gs%ns !ms, ml
- if (wl(i) .gt. anox_crit(i)) then
- qebul(i) = bull(i)*(q_old(i) - ch4_crit(i)) ! cthresh -> ch4_crit(i), according to explicit scheme
- endif
- end do
-
-!Methane bubble flux at the bottom
- febul0 = 0.
- do i = 2, gs%ns-1 !ms+1, ml
-! if(wl(i) .gt. anox_crit(i)) then
-! if(z(i).ge.w) then
-! The original statement is replaced by the grid-consistent one
-! febul = febul + 0.5*(z(i)-z(i-1))*(qebul(i)+qebul(i-1))
- febul0 = febul0 + 0.5*(z(i+1)-z(i-1))*qebul(i)
-! end if
- end do
- febul0 = febul0 + qebul(1)*0.5*(z(2)-z(1))
- febul0 = febul0 + qebul(gs%ns)*0.5*(z(gs%ns)-z(gs%ns-1))
fplant = 0.
- do i = 2, gs%ns-1 !ms+1, ml-1
- fplant = fplant + 0.5*(z(i+1)-z(i-1))*plant(i)*0.5*(qsoil(i) + q_old(i)) * (1 - pox) !* veg
- end do
+ if (soilswitch%par == 1) then
+ do i = 2, gs%ns-1 !ms+1, ml-1
+ fplant = fplant + 0.5*(z(i+1)-z(i-1))*plant(i)*0.5*(qsoil(i) + q_old(i)) * (1 - pox) !* veg
+ end do
+ endif
ifdeep : if (deepestsoil) then
- ! Adding bubble flux from the lowest soil column to horizontally averaged bubble flux
- if (gs%nsoilcols > 1) then
- call BUBBLEFLUXAVER(gs%ix,gs%iy,gs%nsoilcols)
- else
- fbbleflx_ch4_sum(0) = febul0*fbbleflx_ch4(0) *bathymwater(1) %area_int
- fbbleflx_ch4_sum(1:gs%M) = febul0*fbbleflx_ch4(1:gs%M)*bathymwater(1:gs%M)%area_half
- fbbleflx_ch4_sum(gs%M+1) = febul0*fbbleflx_ch4(gs%M+1)*bathymwater(gs%M+1)%area_int
- endif
+ !! Adding bubble flux from the lowest soil column to horizontally averaged bubble flux
+ !if (gs%nsoilcols > 1) then
+ ! call BUBBLEFLUXAVER(gs%ix,gs%iy,gs%nsoilcols)
+ !else
+ ! fbbleflx_ch4_sum(0) = febul0*fbbleflx_ch4(0) *bathymwater(1) %area_int
+ ! fbbleflx_ch4_sum(1:gs%M) = febul0*fbbleflx_ch4(1:gs%M)*bathymwater(1:gs%M)%area_half
+ ! fbbleflx_ch4_sum(gs%M+1) = febul0*fbbleflx_ch4(gs%M+1)*bathymwater(gs%M+1)%area_int
+ !endif
!The second step of splitting-up scheme - bubble dissolution source in water column
if (h1 > 0.) then
@@ -994,9 +1023,12 @@ endif ifdeep
! & diff(2)*(qsoil(2) - qsoil(1))*2./(z(2)-z(1)) + &
! & xx*( - qebul(1) - plant(1)*qsoil(1) + rprod(1) - oxid(1)*qsoil(1))
-
- fdiff1 = - 0.5*diff(gs%ns)*(qsoil(gs%ns) + q_old(gs%ns) - qsoil(gs%ns-1) - q_old(gs%ns-1))/ & ! ml
- & (z(gs%ns) - z(gs%ns-1))
+ if (soilswitch%par == 1) then
+ fdiff1 = - 0.5*diff(gs%ns)*(qsoil(gs%ns) + q_old(gs%ns) - qsoil(gs%ns-1) - q_old(gs%ns-1))/ & ! ml
+ & (z(gs%ns) - z(gs%ns-1))
+ else
+ fdiff1 = 0.
+ endif
! ===============================================================
! total methane emission ftot
@@ -1008,6 +1040,13 @@ endif ifdeep
ftot = fdiff_lake_surf + fbbleflx_ch4_sum(0)/bathymwater(1)%area_int + fplant
endif
+ ! Water methane content
+ do i = 1, gs%M+1
+ q_sum_new = q_sum_new + qwater(i,2)*h1*ddz05(i-1)
+ enddo
+
+ if (soilswitch%par == 1) then
+
q_sum_new = 0.
do i = 2, gs%ns-1 !ms+1, ml-1
q_sum_new = q_sum_new + qsoil(i)*0.5*(z(i+1)-z(i-1))
@@ -1015,11 +1054,6 @@ endif ifdeep
q_sum_new = q_sum_new + qsoil(gs%ns)*0.5*(z(gs%ns)-z(gs%ns-1)) ! ml
q_sum_new = q_sum_new + qsoil(1)*(0.5*(z(2)-z(1))) ! + 0.5*ddz(M)*h1/por(1))
- ! Water methane content
- do i = 1, gs%M+1
- q_sum_new = q_sum_new + qwater(i,2)*h1*ddz05(i-1)
- enddo
-
do i = 2, gs%ns-1 !ms+1, ml-1
rprod_sum = rprod_sum + rprod(i)*dt*0.5*(z(i+1)-z(i-1))
oxid_sum = oxid_sum + oxid(i)*(qsoil(i) + q_old(i))*dt*0.25*(z(i+1)-z(i-1))
@@ -1037,6 +1071,8 @@ endif ifdeep
bull_sum = bull_sum + bull(1)*(q_old(1)-ch4_crit(1))*dt*0.5*(z(2)-z(1))
plant_sum = plant_sum + plant(1)*(qsoil(1) + q_old(1))*dt*0.25*(z(2)-z(1))
+ endif
+
! print*, 'bounds_talik', 0.5*(z(i_talik-1)+z(i_talik)), 0.5*(z(i_talik+1)+z(i_talik))
!if( (q_sum_new - q_sum_old).ne.0. .and. gs%isoilcol == gs%nsoilcols) then
diff --git a/source/model/oxygen_mod.f90 b/source/model/oxygen_mod.f90
index cdb531367f0fa1e6d21aa1520c2fe39755ca6061..2da414866e3a2da13345a83740375e3d0793fa93 100644
--- a/source/model/oxygen_mod.f90
+++ b/source/model/oxygen_mod.f90
@@ -183,11 +183,9 @@ if (ls%h1 > 0) then
! do i = 2, M
! x = x + oxyg(i,2)*(ddz(i-1)+ddz(i))/2.
! enddo
-
endif
-
deallocate (a, b, c, f, y)
END SUBROUTINE OXYGEN