Skip to content
Snippets Groups Projects
Commit 61ca1f2d authored by Victor Stepanenko's avatar Victor Stepanenko
Browse files

Methane simulation in water now works if soilswitch=1 or 0

parent e3d91191
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
source/model/init_var_mod.f90
+ 5
− 5
View file @ 61ca1f2d
......@@ -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
......
source/model/lake.f90
+ 103
− 110
View file @ 61ca1f2d
......@@ -989,6 +989,16 @@ if (soilswitch%par == 1) then
& bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso)
endif
!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
!& fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!& plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!& anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
!& 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
endif
gs%isoilcol = nsoilcols
call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
......@@ -1003,20 +1013,7 @@ if (soilswitch%par == 1) then
& 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
!& fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!& plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!& anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
!& 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
!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,6 +1365,16 @@ if (flag_snow == 1) then
& bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso)
endif
! 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
gs%isoilcol = nsoilcols
call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
......@@ -1382,20 +1389,6 @@ if (flag_snow == 1) then
& 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
!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,6 +1625,16 @@ else
& bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso)
endif
! 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
gs%isoilcol = nsoilcols
call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
......@@ -1646,20 +1649,7 @@ else
& 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
!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,6 +1819,16 @@ if (flag_snow == 1) then
& snowmass_init,a,b,c,d,Temp,phi,fetch,dt)
if (soilswitch%par == 1) then
call BUBBLE_BLOCK
! 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
gs%isoilcol = nsoilcols
call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
......@@ -1843,16 +1843,7 @@ if (flag_snow == 1) then
& 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
l2 = l1 + dhi
h2 = h1 + dhw
......@@ -1910,6 +1901,16 @@ else
if (soilswitch%par == 1) then
call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
call BUBBLE_BLOCK
! 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
gs%isoilcol = nsoilcols
call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
......@@ -1924,16 +1925,6 @@ else
& 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
l2 = l1 + dhi
h2 = h1 + dhw
......@@ -1969,6 +1960,16 @@ if4: IF (layer_case==4) THEN
& snowmass_init,a,b,c,d,Temp,phi,fetch,dt)
if (soilswitch%par == 1) then
call BUBBLE_BLOCK
! 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
gs%isoilcol = nsoilcols
call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
......@@ -1983,16 +1984,7 @@ if4: IF (layer_case==4) THEN
& 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
else
ice = 0; snow = 0; water = 0; deepice = 0
......@@ -2008,6 +2000,16 @@ if4: IF (layer_case==4) THEN
if (soilswitch%par == 1) then
call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
call BUBBLE_BLOCK
! 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
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, &
......@@ -2021,16 +2023,7 @@ if4: IF (layer_case==4) THEN
& 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
! & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
! & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
! & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
! & 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
endif
endif
ENDIF if4
......
source/model/methane_mod.f90
+ 87
− 51
View file @ 61ca1f2d
......@@ -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.
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))
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
......
source/model/oxygen_mod.f90
+ 0
− 2
View file @ 61ca1f2d
......@@ -184,10 +184,8 @@ if (ls%h1 > 0) then
! x = x + oxyg(i,2)*(ddz(i-1)+ddz(i))/2.
! enddo
endif
deallocate (a, b, c, f, y)
END SUBROUTINE OXYGEN
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment