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

Methane equation in water is now can be solved with sediments deactivated,...

Methane equation in water is now can be solved with sediments deactivated, assuming diffusive flux is given at the bottom and ebullition flux is zero
parent d4e7a1b2
No related branches found
No related tags found
No related merge requests found
...@@ -986,6 +986,7 @@ if (soilswitch%par == 1) then ...@@ -986,6 +986,7 @@ if (soilswitch%par == 1) then
& bathymwater,bathymice,bathymdice,bathymsoil, & & bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso) & gsp,o2soilflux, .true.,lso)
endif endif
endif
gs%isoilcol = nsoilcols gs%isoilcol = nsoilcols
call METHANE & call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, & & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
...@@ -1000,20 +1001,7 @@ if (soilswitch%par == 1) then ...@@ -1000,20 +1001,7 @@ if (soilswitch%par == 1) then
& h_talik,tot_ice_meth_bubbles,add_to_winter) & h_talik,tot_ice_meth_bubbles,add_to_winter)
qmethane(1:M+1 ) = qwater(1:M+1,2) qmethane(1:M+1 ) = qwater(1:M+1,2)
qmethane(M+2:M+ns) = qsoil (2:ns,nsoilcols) 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 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, & 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) & o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2)
...@@ -1022,7 +1010,6 @@ call CARBON_DIOXIDE (gs, Tw2, fbbleflx_co2_sum, fbbleflx_co2(0,nsoilcols), lamca ...@@ -1022,7 +1010,6 @@ call CARBON_DIOXIDE (gs, Tw2, fbbleflx_co2_sum, fbbleflx_co2(0,nsoilcols), lamca
call PHOSPHORUS(M,dt,ls,gsp,bathymwater,gas,lsp,lamcarbdi) call PHOSPHORUS(M,dt,ls,gsp,bathymwater,gas,lsp,lamcarbdi)
call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, & call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, &
& qwateroxidtot, qwateroxidML) & qwateroxidtot, qwateroxidML)
!call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), DIC, Tw2) ! two-meth
!Calculation of water density profile at the next timestep !Calculation of water density profile at the next timestep
call DENSITY_W(M,eos%par,lindens%par,Tw2,Sal2,preswat,row2) call DENSITY_W(M,eos%par,lindens%par,Tw2,Sal2,preswat,row2)
...@@ -1367,6 +1354,7 @@ if (flag_snow == 1) then ...@@ -1367,6 +1354,7 @@ if (flag_snow == 1) then
& bathymwater,bathymice,bathymdice,bathymsoil, & & bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso) & gsp,o2soilflux, .true.,lso)
endif endif
endif
gs%isoilcol = nsoilcols gs%isoilcol = nsoilcols
call METHANE & call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, & & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
...@@ -1381,20 +1369,7 @@ if (flag_snow == 1) then ...@@ -1381,20 +1369,7 @@ if (flag_snow == 1) then
& h_talik,tot_ice_meth_bubbles,add_to_winter) & h_talik,tot_ice_meth_bubbles,add_to_winter)
qmethane(1:M+1) = qwater(1:M+1,2) qmethane(1:M+1) = qwater(1:M+1,2)
qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols) 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 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, & 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) & o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2)
...@@ -1403,7 +1378,6 @@ if (flag_snow == 1) then ...@@ -1403,7 +1378,6 @@ if (flag_snow == 1) then
call PHOSPHORUS(M,dt,ls,gsp,bathymwater,gas,lsp,lamcarbdi) call PHOSPHORUS(M,dt,ls,gsp,bathymwater,gas,lsp,lamcarbdi)
call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, & call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, &
& qwateroxidtot, qwateroxidML) & qwateroxidtot, qwateroxidML)
! call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), DIC, Tw2) ! two-meth
!Calculation of water current's velocities and turbulent characteristics !Calculation of water current's velocities and turbulent characteristics
if (Turbpar%par /= 1) then if (Turbpar%par /= 1) then
...@@ -1633,6 +1607,7 @@ else ...@@ -1633,6 +1607,7 @@ else
& bathymwater,bathymice,bathymdice,bathymsoil, & & bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,o2soilflux, .true.,lso) & gsp,o2soilflux, .true.,lso)
endif endif
endif
gs%isoilcol = nsoilcols gs%isoilcol = nsoilcols
call METHANE & call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, & & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
...@@ -1647,20 +1622,7 @@ else ...@@ -1647,20 +1622,7 @@ else
& h_talik,tot_ice_meth_bubbles,add_to_winter) & h_talik,tot_ice_meth_bubbles,add_to_winter)
qmethane(1:M+1) = qwater(1:M+1,2) qmethane(1:M+1) = qwater(1:M+1,2)
qmethane(M+2:M+ns) = qsoil (2:ns,nsoilcols) 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 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, & 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) & o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2)
...@@ -1669,7 +1631,6 @@ else ...@@ -1669,7 +1631,6 @@ else
call PHOSPHORUS(M,dt,ls,gsp,bathymwater,gas,lsp,lamcarbdi) call PHOSPHORUS(M,dt,ls,gsp,bathymwater,gas,lsp,lamcarbdi)
call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, & call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, &
& qwateroxidtot, qwateroxidML) & qwateroxidtot, qwateroxidML)
! call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), DIC, Tw2) ! two-meth
!Calculation of water current's velocities and turbulent characteristics !Calculation of water current's velocities and turbulent characteristics
if (Turbpar%par /= 1) then if (Turbpar%par /= 1) then
...@@ -1830,6 +1791,7 @@ if (flag_snow == 1) then ...@@ -1830,6 +1791,7 @@ if (flag_snow == 1) then
& snowmass_init,a,b,c,d,Temp,phi,fetch,dt) & snowmass_init,a,b,c,d,Temp,phi,fetch,dt)
if (soilswitch%par == 1) then if (soilswitch%par == 1) then
call BUBBLE_BLOCK call BUBBLE_BLOCK
endif
gs%isoilcol = nsoilcols gs%isoilcol = nsoilcols
call METHANE & call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, & & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
...@@ -1844,16 +1806,7 @@ if (flag_snow == 1) then ...@@ -1844,16 +1806,7 @@ if (flag_snow == 1) then
& h_talik,tot_ice_meth_bubbles,add_to_winter) & h_talik,tot_ice_meth_bubbles,add_to_winter)
qmethane(1:M+1) = qwater(1:M+1,2) qmethane(1:M+1) = qwater(1:M+1,2)
qmethane(M+2:M+ns) = qsoil (2:ns,nsoilcols) 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 l2 = l1 + dhi
h2 = h1 + dhw h2 = h1 + dhw
...@@ -1913,6 +1866,7 @@ else ...@@ -1913,6 +1866,7 @@ else
if (saltice%par == 1) call UPDATE_ICE_SALINITY(dt,ls,wst) if (saltice%par == 1) call UPDATE_ICE_SALINITY(dt,ls,wst)
if (soilswitch%par == 1) then if (soilswitch%par == 1) then
call BUBBLE_BLOCK call BUBBLE_BLOCK
endif
gs%isoilcol = nsoilcols gs%isoilcol = nsoilcols
call METHANE & call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, & & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
...@@ -1927,16 +1881,6 @@ else ...@@ -1927,16 +1881,6 @@ else
& h_talik,tot_ice_meth_bubbles,add_to_winter) & h_talik,tot_ice_meth_bubbles,add_to_winter)
qmethane(1:M+1) = qwater(1:M+1,2) qmethane(1:M+1) = qwater(1:M+1,2)
qmethane(M+2:M+ns) = qsoil (2:ns,nsoilcols) 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 l2 = l1 + dhi
h2 = h1 + dhw h2 = h1 + dhw
...@@ -1972,6 +1916,7 @@ if4: IF (layer_case==4) THEN ...@@ -1972,6 +1916,7 @@ if4: IF (layer_case==4) THEN
& snowmass_init,a,b,c,d,Temp,phi,fetch,dt) & snowmass_init,a,b,c,d,Temp,phi,fetch,dt)
if (soilswitch%par == 1) then if (soilswitch%par == 1) then
call BUBBLE_BLOCK call BUBBLE_BLOCK
endif
gs%isoilcol = nsoilcols gs%isoilcol = nsoilcols
call METHANE & call METHANE &
& (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, & & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
...@@ -1986,16 +1931,6 @@ if4: IF (layer_case==4) THEN ...@@ -1986,16 +1931,6 @@ if4: IF (layer_case==4) THEN
& h_talik,tot_ice_meth_bubbles,add_to_winter) & h_talik,tot_ice_meth_bubbles,add_to_winter)
qmethane(1:M+1) = qwater(1:M+1,2) qmethane(1:M+1) = qwater(1:M+1,2)
qmethane(M+2:M+ns) = qsoil (2:ns,nsoilcols) 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 else
ice = 0; snow = 0; water = 0; deepice = 0 ice = 0; snow = 0; water = 0; deepice = 0
...@@ -2011,6 +1946,7 @@ if4: IF (layer_case==4) THEN ...@@ -2011,6 +1946,7 @@ if4: IF (layer_case==4) THEN
if (soilswitch%par == 1) then if (soilswitch%par == 1) then
call SOILFORLAKE(dt,a,b,c,d,nsoilcols) call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
call BUBBLE_BLOCK call BUBBLE_BLOCK
endif
gs%isoilcol = nsoilcols gs%isoilcol = nsoilcols
call METHANE(gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, & 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,Sals1, & & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals1, &
...@@ -2024,16 +1960,6 @@ if4: IF (layer_case==4) THEN ...@@ -2024,16 +1960,6 @@ if4: IF (layer_case==4) THEN
& h_talik,tot_ice_meth_bubbles,add_to_winter) & h_talik,tot_ice_meth_bubbles,add_to_winter)
qmethane(1:M+1) = qwater(1:M+1,2) qmethane(1:M+1) = qwater(1:M+1,2)
qmethane(M+2:M+ns) = qsoil (2:ns,nsoilcols) 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
ENDIF if4 ENDIF if4
......
SUBROUTINE METHANE2 &
& (pressure,wind10,zsoil,Tsoil,wl,wi,wa,rootss,rhogr,por,veg,qsoil,TgrAnn, &
& ddz, Twater, lammeth, qwater, h1, l1, ls1, dhw, dhw0, &
& fplant, febul, fdiff, ftot, fdiff_lake_surf, &
& plant_sum,bull_sum,oxid_sum,rprod_sum, &
& anox,M,ns,dt,eps_surf,rprod_total_oldC,rprod_total_newC, &
& h_talik,tot_ice_meth_bubbles)
use LAKE_DATATYPES, only : ireals, iintegers
use DRIVING_PARAMS !, only : &
! & tricemethhydr
use PHYS_CONSTANTS, only : &
& row0, g, row0, roa0, &
& row0_d_roi, row0_d_roa0, &
& roa0_d_roi, roa0_d_row0
use METH_OXYG_CONSTANTS
use PHYS_PARAMETERS
use NUMERIC_PARAMS
use PHYS_FUNC, only : &
& WL_MAX, &
& HENRY_CONST, &
& DIFF_WATER_METHANE, &
& DIFF_AIR_METHANE, &
& GAS_WATATM_FLUX
use NUMERICS, only : &
PROGONKA
implicit none
!* ===============================================================
!* prod + oxid
!* solution 1-D nonlineal equation
!* describing production, oxidation
!* and transport of methane from wetlands
!* ===============================================================
!* input parameters:
!* =====================
!* z(ml) : height of level in soil from soil surface (cm)
!* Tsoil(ml) : the soil temperature (C degrees)
!*
!* output parameters:
!* =====================
!* fplant : the ch4 flux due to plant-mediated transport
!* febul : the ch4 flux due to ebullition
!* fdiff : the diffusive ch4 flux at soil/water=atm boundary
!* ftot : the total ch4 emission
!* =============================================================
!* parameters for production
!* =============================
!* rnpp ! NPP (gC/(m**2*mo))
!* rnppmax ! the maximum value of the NPP
!* rnroot ! the rooting depth (cm)
!* q100 ! value lying within the range of observed
!* ! q100 values ranging from 1.7 to 16
!* r0 ! the constant rate factor (muM/s)
!* =============================================================
!* parameters for oxidation
!* =============================
!* vmax ! Michaelis-Menten (muM/s)
!* rkm ! coefficients (muM)
!* q10 ! the observed value for oxidation
!* TgrAnn(ml) ! the annual mean soil temperature (deg C)
!* =============================================================
!* parameters for ebullition
!* =============================
!* rke ! a rate constant of the unit 1/s
!* cmin ! the concentration at which bubble
!* ! formation occurs (muM)
!* punveg ! the percentage of unvegetated, bare soil
!* cthresh ! the threshold concentration
!* ! for bubble formation
!* =============================================================
!* parameters for plant-mediated transport
!* ===========================================
!* rkp ! a rate constant of the unit 0.01/s
!* tveg ! a factor describing the quality of
!* ! plant-mediated transport at a site
!* pox ! a certain fraction of methane oxidized in
!* ! the oxic zone around the roots of plants
!* rlmin ! the parameter for the function of fgrow(t)
!* rl ! the parameter for the function of fgrow(t)
!* rlmax ! rlmax=rlmin+rl
!* ===============================================================
real(kind=ireals), intent(in) :: pressure ! Atmospheric pressure, Pa
real(kind=ireals), intent(in) :: wind10 ! Wind speed at 10 m above the surface
real(kind=ireals), intent(in) :: zsoil(1:ns) ! numerical levels in soil, meters
real(kind=ireals), intent(in) :: Tsoil(1:ns) ! Celsius
real(kind=ireals), intent(in) :: wl(1:ns) ! liquid water content in a soil, kg/kg
real(kind=ireals), intent(in) :: wi(1:ns) ! ice content in a soil, kg/kg
real(kind=ireals), intent(in) :: wa(1:ns) ! air content in a soil, kg/kg
! COMMON /VEINIT/
real(kind=ireals), intent(in) :: rootss(1:ns) ! meters
! COMMON /SOILDAT/
real(kind=ireals), intent(in) :: rhogr(1:ns) ! density of dry soil, kg/m**3
real(kind=ireals), intent(in) :: por(1:ns) ! soil porosity, n/d
real(kind=ireals), intent(in) :: ddz(1:M) ! the thickness of dzeta-layers in water, n/d
real(kind=ireals), intent(in) :: Twater(1:M+1) ! the water temperature, degrees Celsius
real(kind=ireals), intent(in) :: lammeth(1:M) ! the heat turbulent diffusivity in water, m**2/s
real(kind=ireals), intent(in) :: h1, l1, ls1 ! the thickness (depth) of the water column, m
! ice and bottom ice
real(kind=ireals), intent(in) :: dhw, dhw0 ! the time increment of water column thickness
! and those at its top, m
real(kind=ireals), intent(in) :: dt ! Timestep, sec
real(kind=ireals), intent(in) :: eps_surf ! TKE dissipation rate, m**2/s**3
integer(kind=iintegers), intent(in) :: M, ns ! Number of levels in water and soil
real(kind=ireals), intent(inout) :: veg
! common/ch4/
real(kind=ireals), intent(inout) :: qsoil(1:ns,1:2) ! Methane concentration in soil
real(kind=ireals), intent(inout) :: qwater(1:M+1,1:2,1:2) ! Methane concentration in water
real(kind=ireals), intent(inout) :: TgrAnn(1:ns)
real(kind=ireals), intent(inout) :: fplant, ftot, fdiff_lake_surf(1:2)
real(kind=ireals), intent(inout) :: febul(1:2), fdiff(1:2)
real(kind=ireals), intent(inout) :: plant_sum, bull_sum,oxid_sum,rprod_sum
real(kind=ireals), intent(inout) :: anox
real(kind=ireals), intent(inout) :: tot_ice_meth_bubbles(1:2)
real(kind=ireals), intent(out) :: rprod_total_oldC
real(kind=ireals), intent(out) :: rprod_total_newC
real(kind=ireals), intent(out) :: h_talik ! talik depth, m
! Local variables
integer(kind=iintegers), parameter :: water_methane_indic = 8
real(kind=ireals), parameter :: C_depth_age = 0.25 ! from 0.25 to 0.30, according to West and Plug, 2007
real(kind=ireals), parameter :: C_talik_age = 0.5 ! from 0.5 to 0.7, according to West and Plug, 2007
real(kind=ireals), parameter :: Tmelt_pnt = 0. ! melting point temperature, degrees Celsius
real(kind=ireals), save :: fcoarse = 0.5 !the relative volume of the coarse pores
real(kind=ireals), save :: scf = 16.04/(1./864.)! the scale factor mgCH4/m**2/day
! scf = 12./(24.*1./864)! the scale factor mgC/m**2/h
real(kind=ireals) :: diff(ns), forg(ns)
real(kind=ireals) :: qplant(ns),qebul(1:ns,1:2),q_old(ns),bull(1:ns,1:2)
real(kind=ireals) :: ch4_crit(1:ns,1:2)
real(kind=ireals) :: oxid(ns),dz_soil(ns),z2(ns)
real(kind=ireals) :: fcch4(ns,1:2)
real(kind=ireals) :: bp(ns), ap(ns)
real(kind=ireals) :: rprod(ns), rprod_new(ns), rprod_old(ns), plant(ns)
real(kind=ireals) :: n2(ns)
real(kind=ireals) :: water_porvolrat(1:ns), air_porvolrat(1:ns)
real(kind=ireals), allocatable :: a(:), b(:), c(:), f(:), y(:)
! real(kind=ireals), allocatable, save :: mean_rprod_oldc0(:)
real(kind=ireals) :: rnroot
real(kind=ireals) :: froot(ns), anox_crit(ns)
real(kind=ireals) :: pow_oldorg
real(kind=ireals) :: lake_age
real(kind=ireals) :: tmp_veg, t_grow, t_mature, fnpp
real(kind=ireals) :: fin, ft, ft1, t50, fgrow, q_sum_old, fdiff1, q_sum_new
real(kind=ireals) :: balance
real(kind=ireals) :: Flux_atm(1:2)
real(kind=ireals) :: ratio_c13_c14
real(kind=ireals) :: xx, yy ! working variables
real(kind=ireals), allocatable :: z(:)
real(kind=ireals), allocatable :: roots(:)
real(kind=ireals), external :: DZETA
integer(kind=iintegers) :: i, j ! Loop indices
integer(kind=iintegers) :: i_talik
real(kind=ireals), external :: MEAN_RPROD_OLDC, MEAN_RPROD_OLDC2
allocate (z(1:ns))
allocate (roots(1:ns))
allocate (a(1:vector_length),b(1:vector_length),c(1:vector_length), &
& f(1:vector_length),y(1:vector_length))
z(:) = zsoil(:)
roots(:) = rootss(:)
! TgrAnn(ns) = Tsoil(ns) !ms
TgrAnn(ns) = 0. ! exludes the factor of annual mean temperature
! Calculation of pore volume ratios for liquid water and air
! They are needed for diffusivity calculations
do i = 1, ns
water_porvolrat(i) = wl(i) / &
& (por(i)*(row0/rhogr(i) + wl(i) + row0_d_roi*wi(i) + row0_d_roa0*wa(i)) )
air_porvolrat(i) = wa(i) / &
& (por(i)*(roa0/rhogr(i) + wa(i) + roa0_d_roi*wi(i) + roa0_d_row0*wl(i)) )
enddo
! Assuming exponential dependence of nitrogen concentration on depth
n2(1) = 2.*n2_atm0*(1.-exp(-n2_exp_decay*0.5*z(2)) )/(z(2)*n2_exp_decay)
n2(2:ns) = 0. ! nitrogen concentration in the layers below the surface is
! set to zero, since it decays rapidly with depth (Bazhin, 2001)
! Calculation of the talik depth and of the power in production term due to
! old talik organics decomposition
i_talik = 0
do i = 1, ns-1
if (Tsoil(i) > Tmelt_pnt .and. Tsoil(i+1) <= Tmelt_pnt) then
h_talik = z(i) + (z(i+1) - z(i)) * &
& (Tmelt_pnt - Tsoil(i))/(Tsoil(i+1) - Tsoil(i))
if (h_talik < 0.5*(z(i) + z(i+1))) then
i_talik = i
else
i_talik = i + 1
endif
exit
endif
enddo
! if (.not.allocated(mean_rprod_oldc0)) then
! allocate (mean_rprod_oldc0(1:ns))
! do i = 2, ns-1
! mean_rprod_oldc0(i) = &
! & MEAN_RPROD_OLDC(0.5*(z(i-1) + z(i)), 0.5*(z(i+1) + z(i)), &
! & h_talik, C_talik_age)
! enddo
! mean_rprod_oldc0(1) = &
! & MEAN_RPROD_OLDC(z(1), 0.5*(z(2) + z(1)), &
! & h_talik, C_talik_age)
! mean_rprod_oldc0(ns) = &
! & MEAN_RPROD_OLDC(0.5*(z(ns-1) + z(ns)), z(ns), &
! & h_talik, C_talik_age)
! endif
! Calculation of approximate THERMOKARST lake age, according to results by West and Plug, 2007
! lake_age = (h1/C_depth_age)**2 ! years
! lake_age = (h_talik/C_talik_age)**2 ! years
! pow_oldorg = 2.*alpha_old_org*lake_age
do i = 1, ns
xx = HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, Henry_temp_ref, Tsoil(i)+273.15)
yy = HENRY_CONST(Henry_const0_n2, Henry_temp_dep_n2, Henry_temp_ref, Tsoil(i)+273.15)
ch4_crit(i,1) = qsoil(i,1)*rel_conc_ebul_crit*por(i)*xx*(pressure + row0*g*h1 - n2(i)/yy) / &
& (qsoil(i,1) + qsoil(i,2))
ch4_crit(i,2) = ch4_crit(i,1)*qsoil(i,2)/qsoil(i,1)
! pressure -> 1.d+5
enddo
do i = 1, ns !ms+1, ml
! anox_crit(i) = (por(i)*rhow/rhogr(i)-wi(i)*rhow/rhoi)*0.9
anox_crit(i) = WL_MAX(por(i),rhogr(i),wi(i),tricemethhydr%par)*0.9
end do
! Calculation of rooting depth
rnroot = 0.
do i = 1, ns !ms, ml
if (roots(i).gt.0.) rnroot = z(i)
end do
!* ============================================================
!* tgrow is the temperature at which plants start to grow
!* ============================================================
if (Tsoil(2) .lt. 5.) then ! 5 deg. Celsius
t_grow = 2. ! deg. Celsius
else
t_grow = 7. ! deg. Celsius
end if
t_mature = t_grow + 10. ! the temperature at which plants reach maturity, deg. Celsius
!* =============================================================
!* parameters for diffusion
!* =============================
anox = 0.
do i = ns-1, 2, -1 !ml-1, ms+1, -1
if (z(i) .le. 1.) then ! 100 cm
if (wl(i) .ge. anox_crit(i)) then
anox = anox + (z(i+1)-z(i-1))*0.5
end if
end if
end do
if (rnroot > 0.) then
do i = 1, ns !ms, ml
if(z(i) .le. rnroot) then
froot(i) = 2.*(rnroot-z(i))/rnroot
else
froot(i) = 0.
end if
end do
else
froot(:) = 0.
endif
do i = 1, ns
if (z(i) .gt. rnroot) then ! for vegetated
forg(i) = veg * exp((rnroot - z(i))*10.) ! soil
else
! if (z(i).le.rnroot) then
forg(i) = veg * 1.
endif
end do
do i = 1, ns ! for unvegetated
forg(i) = forg(i) + (1.-veg) * forg0*exp(-z(i)*lambda_new_org) ! soils
! forg(i) = forg(i) + (1.-veg) * forg0 ! the option for thermokarst lakes
! for which it is assumed the organic
! matter is homogeneously distributed
! in soil
end do
do i = 2, ns
! if (wl(i) .lt. anox_crit(i)) then
! diff(i) = diff_unsat*tortuosity_coef*fcoarse ! in the unsaturated
!
! else ! soil layers
! diff(i) = diff_water*tortuosity_coef*fcoarse ! in the water saturated
! diff(i) = DIFF_WATER_METHANE(Tsoil(i)) * &
! & tortuosity_coef*fcoarse
! end if ! soil layers
diff(i) = 0.25*( por(i) + por(i-1) )*tortuosity_coef * &
& ( (water_porvolrat(i) + water_porvolrat(i-1)) * &
& DIFF_WATER_METHANE( 0.5*(Tsoil(i) + Tsoil(i-1)) ) + &
& (air_porvolrat(i) + air_porvolrat(i-1)) * &
& DIFF_AIR_METHANE( 0.5*(Tsoil(i) + Tsoil(i-1)) ) )
end do
fnpp = rnpp !function of the variation of NPP with t
fin = 1. + fnpp/rnppmax !variation of substrate availability with t
rprod_total_newC = 0.
rprod_total_oldC = 0.
do i = 1, ns
if (wl(i) .ge. anox_crit(i)) then
ft1 = 1. ! Note, methane generation through old organics decomposition
! is zero below talik automatically, see the code below
if (Tsoil(i) .gt. 0) then
if (i /= i_talik) then
ft = 1.
else
ft = (h_talik - 0.5*(z(i-1) + z(i)) ) / &
& (0.5*(z(i+1) - z(i-1)))
endif
else
if (i /= i_talik) then
ft = 0.
else
ft = (h_talik - 0.5*(z(i-1) + z(i)) ) / &
& (0.5*(z(i+1) - z(i-1)))
endif
endif
else
ft = 0.
if (i /= ns .and. i /= 1) then
if ((Tsoil(i) > 0. .and. Tsoil(i+1) < 0.) .or. &
& (Tsoil(i) < 0. .and. Tsoil(i-1) > 0.) ) then
! Regularization since there is always minimum of soil moisture content
! near the talik depth (due to diffusion from water saturated layers above
! to frozen layers below)
ft1 = 1. - (anox_crit(i) - wl(i))/anox_crit(i)
! ft1 = 1.
! ft = ft1
! if (i == i_talik) &
! & ft = ft * (h_talik - 0.5*(z(i-1) + z(i)) ) / &
! & (0.5*(z(i+1) - z(i-1)))
else
ft1 = 0.
endif
else
ft1 = 0.
endif
endif
! xx = r0_oldorg*(min(zsoil(i)/h_talik,1.d0))**pow_oldorg * &
! & q100**((Tsoil(i)-TgrAnn(ns))*0.1)*ft
! xx = r0_oldorg*exp(-alpha_old_org*C_talik_age**(-2) * &
! & (h_talik**2-(min(z(i),h_talik))**2))*ft * &
! & q100**((Tsoil(i)-TgrAnn(ns))*0.1)
ratio_c13_c14 = min(z(i)/h_talik, 1.d0)
! Note that for the performance optimization MEAN_RPROD_OLDCs may
! be computed only once - at the first timestep
! Note that calculating MEAN_RPROD_OLDC the temperature step
! function weighting is already implemented, so ft multiplication
! is not needed
if (i == 1) then
! xx = MEAN_RPROD_OLDC(0.d0, 0.5*(z(2) + z(1)), h_talik, C_talik_age)
xx = MEAN_RPROD_OLDC2(0.d0, 0.5*(z(2) + z(1)), h_talik, C_talik_age)
! xx = exp(-alpha_old_org*C_talik_age**(-2)*h_talik**2) * &
! & mean_rprod_oldc0(1)
rprod_old(i) = xx*q100**((Tsoil(1)-TgrAnn(ns))*0.1)*ft1*ratio_c13_c14
rprod_new(i) = xx*q100**((Tsoil(1)-TgrAnn(ns))*0.1)*ft1*(1. - ratio_c13_c14)
elseif (i == ns) then
! xx = &
! & ( r0_oldorg * &
! & exp(-alpha_old_org*C_talik_age**(-2) * &
! & (h_talik**2-min(0.5*(z(ns)+z(ns-1)),h_talik)**2 ) ) + &
! & r0_oldorg * &
! & exp(-alpha_old_org*C_talik_age**(-2) * &
! & (h_talik**2-(min(z(ns),h_talik))**2 ) ) )* 0.5
! xx = MEAN_RPROD_OLDC(0.5*(z(ns) + z(ns-1)), z(ns), h_talik, C_talik_age)
xx = MEAN_RPROD_OLDC2(0.5*(z(ns) + z(ns-1)), z(ns), h_talik, C_talik_age)
! xx = exp(-alpha_old_org*C_talik_age**(-2)*h_talik**2) * &
! & mean_rprod_oldc0(ns)
rprod_old(i) = xx*q100**((Tsoil(ns)-TgrAnn(ns))*0.1)*ft1*ratio_c13_c14
rprod_new(i) = xx*q100**((Tsoil(ns)-TgrAnn(ns))*0.1)*ft1*(1. - ratio_c13_c14)
else
! xx = &
! & ( r0_oldorg * &
! & exp(-alpha_old_org*C_talik_age**(-2) * &
! & (h_talik**2-(min(0.5*(z(i)+z(i-1)),h_talik))**2)) + &
! & r0_oldorg * &
! & exp(-alpha_old_org*C_talik_age**(-2) * &
! & (h_talik**2-(min(z(i),h_talik))**2)) ) * 0.5 * &
! & (z(i)-z(i-1))/(z(i+1)-z(i-1)) + &
! ( r0_oldorg * &
! & exp(-alpha_old_org*C_talik_age**(-2) * &
! & (h_talik**2-(min(0.5*(z(i+1)+z(i)),h_talik))**2)) + &
! & r0_oldorg * &
! & exp(-alpha_old_org*C_talik_age**(-2) * &
! & (h_talik**2-(min(z(i),h_talik))**2)) ) * 0.5 * &
! & (z(i+1)-z(i))/(z(i+1)-z(i-1))
! xx = MEAN_RPROD_OLDC(0.5*(z(i) + z(i-1)), 0.5*(z(i+1) + z(i)), h_talik, C_talik_age)
xx = MEAN_RPROD_OLDC2(0.5*(z(i) + z(i-1)), 0.5*(z(i+1) + z(i)), h_talik, C_talik_age)
! xx = exp(-alpha_old_org*C_talik_age**(-2)*h_talik**2) * &
! & mean_rprod_oldc0(i)
rprod_old(i) = xx*q100**((Tsoil(i)-TgrAnn(ns))*0.1)*ft1*ratio_c13_c14
rprod_new(i) = xx*q100**((Tsoil(i)-TgrAnn(ns))*0.1)*ft1*(1. - ratio_c13_c14)
endif
rprod_new(i) = rprod_new(i) + &
& r0_oliglake*forg(i)*fin * &
& q100**((Tsoil(i)-TgrAnn(ns))*0.1)*ft
if (i == 1) then
rprod_total_oldC = rprod_total_oldC + rprod_old(1)*0.5*(z(2)-z(1))
rprod_total_newC = rprod_total_newC + rprod_new(1)*0.5*(z(2)-z(1))
elseif (i == ns) then
rprod_total_oldC = rprod_total_oldC + rprod_old(ns)*0.5*(z(ns)-z(ns-1))
rprod_total_newC = rprod_total_newC + rprod_new(ns)*0.5*(z(ns)-z(ns-1))
else
rprod_total_oldC = rprod_total_oldC + rprod_old(i)*0.5*(z(i+1)-z(i-1))
rprod_total_newC = rprod_total_newC + rprod_new(i)*0.5*(z(i+1)-z(i-1))
endif
! rprod(i) = rprod_new(i) + rprod_old(i)
! if(qsoil(i).ge.cthresh) then
do j = 1, 2
if (qsoil(i,j) .ge. ch4_crit(i,j)) then
fcch4(i,j) = 1.
else
fcch4(i,j) = 0.
end if
enddo
end do
!* ===============================================================
!* the function fgrow describes the growing state of the plants
!* ===============================================================
!* t50 - the soil temperature at 50cm depth below soil surface
do i = 1, ns !ms, ml
if (z(i) .le. 0.5) then
t50 = Tsoil(i)
else
exit
endif
end do
if (t50.lt.t_grow) fgrow = rlmin
if (t50 .ge. t_grow .and. t50 .le. t_mature) then
fgrow = rlmin + rl*(1.-((t_mature-t50)/(t_mature-t_grow))**2)
endif
if (t50.gt.t_mature) fgrow = rlmax
! =============================================================
do i = 1, ns !ms, ml
!* oxid(i) - coefficient for oxidation
!* methane oxidation occurs only in the unsaturated soil
!* layers
if (wl(i).lt.anox_crit(i)) then
! if(z(i).lt.w) then
! oxid(i) = vmax*q10**((Tsoil(i) - TgrAnn(ns))*0.1)/(rkm + qsoil(i)) !ml
! Currently oxidation is turned off in soil layers.
! It should be taken into account when oxygen concentration will
! be somehow extimated (calculated).
oxid(i) = 0.
else
oxid(i) = 0.
end if
!* plant(i) - coefficient for plant-mediated transport
! plant(i) = rkp*tveg*froot(i)*fgrow*(1.-0.5)*veg
! plant(i) = rkp*tveg*froot(i)*fgrow*veg
plant(i) = 0.
!*
!* bull(i) - coefficient for ebullition
! if(z(i).ge.w) then
do j = 1, 2
if (wl(i) .ge. anox_crit(i)) then
bull(i,j) = rke*fcch4(i,j)
else
bull(i,j) = 0.
end if
enddo
end do
! write (*,*) bull(:,2)
! write (*,*) rprod_old(:)
! write (*,*) ch4_crit(:,2)
! oxid_sum = 0.
! rprod_sum = 0.
! bull_sum = 0.
! plant_sum = 0.
! q_sum_old = 0.
! do i = 2, ns-1 !ms+1, ml-1
! q_sum_old = q_sum_old + qsoil(i)*0.5*(z(i+1)-z(i-1))
! end do
!
! q_sum_old = q_sum_old + qsoil(ns)*0.5*(z(ns)-z(ns-1))
! q_sum_old = q_sum_old + qsoil(1)*(0.5*(z(2)-z(1)) + 0.5*ddz(M)*h1)
! do i = 1, ns !1, ml
! q_old(i) = qsoil(i,j)
! end do
!* =============================================================
!* set up the tridiagonal system of equations
!* =============================================================
meth2 : do j = 1, 2
if (j == 1) rprod(:) = rprod_old(:)
if (j == 2) rprod(:) = rprod_new(:)
if (h1 > 0) then
if (l1 == 0) then
! c(1) = 1.
! b(1) = 0.
! f(1) = ch4_pres_atm0 * &
! & HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, Henry_temp_ref, Twater(1)+273.15)
! ch4_atm0 ! atmospheric concentration
xx = HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, Henry_temp_ref, Twater(1)+273.15)
Flux_atm(j) = GAS_WATATM_FLUX &
& (Twater(1),wind10,qwater(1,1,j),ch4_pres_atm0,xx,water_methane_indic,eps_surf)
else
Flux_atm(j) = 0.
endif
c(1) = - lammeth(1)/(h1*ddz(1)) + dhw0/(2.d0*dt) - &
& ddz(1)*h1/(2.d0*dt)
b(1) = - lammeth(1)/(h1*ddz(1)) + dhw0/(2.d0*dt)
f(1) = - ddz(1)*h1*qwater(1,1,j)/(2.d0*dt) + Flux_atm(j)
call DIFF_COEF(a,b,c,f,2,M,2,M,water_methane_indic,dt)
if (ls1 > 0) then
! Case water, deepice and soil; upper ice and snow are allowed
! Methane diffusion in water
c(M+1) = - lammeth(M)/(ddz(M)*h1)-ddz(M)*h1/(2.d0*dt) + & ! bug corrected: - (dhw-dhw0)/(2.d0*dt)
& (dhw-dhw0)/(2.d0*dt)
a(M+1) = - lammeth(M)/(ddz(M)*h1)+(dhw-dhw0)/(2.d0*dt)
f(M+1) = - qwater(M+1,1,j)*ddz(M)*h1/(2.d0*dt) ! + qflux1
call PROGONKA (a,b,c,f,y,1,M+1)
y(1:M+1) = max(y(1:M+1),0.d0)
qwater(1:M+1,2,j) = y(1:M+1)
! Methane evolution in soil under bottom ice
b(1) = 2.*dt*diff(2)/((z(2)-z(1))*(z(2)-z(1)))
c(1) = 1. + b(1) + dt*(bull(1,j) + plant(1) + oxid(1))
! No diffusive flux is assumed at the ice-soil interface
f(1) = qsoil(1,j) + dt*(rprod(1) + bull(1,j)*ch4_crit(1,j)) ! cthresh -> ch4_crit(1)
do i = 2, ns-1 !ms+1, ml-1
ap(i) = 2.*dt/((z(i+1)-z(i-1))*(z(i+1)-z(i)))
bp(i) = 2.*dt/((z(i+1)-z(i-1))*(z(i)-z(i-1)))
a(i) = bp(i) * diff(i)
c(i) = 1. + ap(i)*diff(i+1) + bp(i)*diff(i) + &
& dt * (bull(i,j)+plant(i)+oxid(i))
f(i) = qsoil(i,j) + dt*rprod(i) + dt*bull(i,j)*ch4_crit(i,j) ! cthresh -> ch4_crit(i)
b(i) = ap(i)*diff(i+1)
end do
bp(ns) = 2.*dt/( (z(ns)-z(ns-1))*(z(ns)-z(ns-1)) )
a(ns) = bp(ns) * diff(ns) ! diff(ns-1)
c(ns) = 1. + a(ns) + & ! bp(ns)*diff(ns)
& dt * (bull(ns,j) + plant(ns) + oxid(ns))
f(ns) = qsoil(ns,j) + dt*(rprod(ns) + bull(ns,j)*ch4_crit(ns,j) ) ! cthresh -> ch4_crit(ns)
call PROGONKA(a, b, c, f, y, 1, ns)
y(1:ns) = max(y(1:ns),0.d0)
qsoil(1:ns,j) = y(1:ns)
else
! Case water, soil; upper ice and snow are allowed
!---------------------WATER-SOIL INTERFACE------------------
a(M+1)= lammeth(M)/(ddz(M)*h1) - (dhw-dhw0)/(2.*dt)
b(M+1)= diff(2)/( (z(2)-z(1)) )
c(M+1) = a(M+1) + b(M+1) + (z(2)-z(1))/(2.*dt) + &
& ddz(M)*h1/(2.*dt) + 0.5*(z(2)-z(1))*(bull(1,j) + plant(1) + oxid(1))
f(M+1)= qwater(M+1,1,j)*((z(2)-z(1))/(2.*dt) + &
& ddz(M)*h1/(2.*dt)) + 0.5*(z(2)-z(1))*(rprod(1) + bull(1,j)*ch4_crit(1,j)) ! cthresh -> ch4_crit(1)
!-----------------------------------------------------------
do i = 2, ns-1 !ms+1, ml-1
ap(i) = 2.*dt/((z(i+1)-z(i-1))*(z(i+1)-z(i)))
bp(i) = 2.*dt/((z(i+1)-z(i-1))*(z(i)-z(i-1)))
a(i+M) = bp(i) * diff(i)
c(i+M) = 1. + ap(i)*diff(i+1) + bp(i)*diff(i) + &
& dt * (bull(i,j)+plant(i)+oxid(i))
f(i+M) = qsoil(i,j) + dt*rprod(i) + dt*bull(i,j)*ch4_crit(i,j) ! cthresh -> ch4_crit(i)
b(i+M) = ap(i)*diff(i+1)
end do
bp(ns) = 2.*dt/( (z(ns)-z(ns-1))*(z(ns)-z(ns-1)) )
a(M+ns) = bp(ns) * diff(ns) ! diff(ns-1)
c(M+ns) = 1. + a(M+ns) + & ! bp(ns)*diff(ns)
& dt * (bull(ns,j) + plant(ns) + oxid(ns))
f(M+ns) = qsoil(ns,j) + dt*(rprod(ns) + bull(ns,j)*ch4_crit(ns,j)) ! cthresh -> ch4_crit(ns)
call PROGONKA (a,b,c,f,y,1,M+ns)
y(1:M+ns) = max(y(1:M+ns),0.d0)
qwater(1:M+1,2,j) = y(1:M+1)
qsoil(1:ns,j) = y(M+1:M+ns)
endif
else
! Case of the soil and the ice above
! Methane evolution in soil under the ice
b(1) = 2.*dt*diff(2)/((z(2)-z(1))*(z(2)-z(1)))
c(1) = 1. + b(1) + dt*(bull(1,j) + plant(1) + oxid(1))
! No diffusive flux is assumed at the ice-soil interface
f(1) = qsoil(1,j) + dt*(rprod(1) + bull(1,j)*ch4_crit(1,j)) ! cthresh -> ch4_crit(1)
do i = 2, ns-1 !ms+1, ml-1
ap(i) = 2.*dt/((z(i+1)-z(i-1))*(z(i+1)-z(i)))
bp(i) = 2.*dt/((z(i+1)-z(i-1))*(z(i)-z(i-1)))
a(i) = bp(i) * diff(i)
c(i) = 1. + ap(i)*diff(i+1) + bp(i)*diff(i) + &
& dt * (bull(i,j)+plant(i)+oxid(i))
f(i) = qsoil(i,j) + dt*rprod(i) + dt*bull(i,j)*ch4_crit(i,j) ! cthresh -> ch4_crit(i)
b(i) = ap(i)*diff(i+1)
end do
bp(ns) = 2.*dt/( (z(ns)-z(ns-1))*(z(ns)-z(ns-1)) )
a(ns) = bp(ns) * diff(ns) ! diff(ns-1)
c(ns) = 1. + a(ns) + & ! bp(ns)*diff(ns)
& dt * (bull(ns,j) + plant(ns) + oxid(ns))
f(ns) = qsoil(ns,j) + dt*rprod(ns) + dt*bull(ns,j)*ch4_crit(ns,j) ! cthresh -> ch4_crit(ns)
call PROGONKA(a, b, c, f, y, 1, ns)
y(1:ns) = max(y(1:ns),0.d0)
qsoil(1:ns,j) = y(1:ns)
endif
enddo meth2
!* ============================================================
!* the methane flux due to plant-mediated transport
!* =================fplant(t)==================================
! fplant = 0.
! do i = 2, ns-1 !ms+1, ml-1
! fplant = fplant + 0.5*(z(i+1)-z(i-1))*plant(i)*qsoil(i) * (1 - pox) !* veg
! end do
!* =============================================================
!* the methane flux due to ebullition
!* =================febul(t)====================================
do j = 1, 2
do i = 1, ns !ms, ml
qebul(i,j) = bull(i,j)*(qsoil(i,j) - ch4_crit(i,j)) ! cthresh -> ch4_crit(i)
end do
enddo
febul(:) = 0.
do j = 1, 2
do i = 2, ns-1
if (wl(i) .gt. anox_crit(i)) then
febul(j) = febul(j) + 0.5*(z(i+1)-z(i-1))*qebul(i,j)
end if
enddo
febul(j) = febul(j) + qebul(1,j)*0.5*(z(2)-z(1))
febul(j) = febul(j) + qebul(ns,j)*0.5*(z(ns)-z(ns-1))
end do
if (l1 /= 0.) then
do j = 1, 2
! A portion of methane ebullition flux (meth_ebul_ice_intercept)
! is intercepted by the ice cover
febul(j) = febul(j)*(1. - meth_ebul_ice_intercept)
tot_ice_meth_bubbles(j) = tot_ice_meth_bubbles(j) + febul(j)*meth_ebul_ice_intercept*dt
enddo
endif
if (l1 == 0.) then
do j = 1, 2
if (tot_ice_meth_bubbles(j) /= 0.) then
! When ice cover disappears, all ice-trapped bubbles are
! released to the atmosphere instantaneously
febul(j) = febul(j) + tot_ice_meth_bubbles(j)/dt
tot_ice_meth_bubbles(j) = 0.
endif
enddo
endif
! ===============================================================
!* the diffusive flux at the soil/water boundary, downwards
!* =================fdiff(t, dzeta=1)==================================
! fdiff = - (diff(1) + diff(2))*0.5*(qsoil(1)-qsoil(2))/ & ! ms
! & (z(2)-z(1))
! fdiff = 0.5*h1*ddz(M)/dt*(qwater(M+1,2)-qwater(M+1,1)) - &
! & lammeth(M)*(qwater(M+1,2)-qwater(M,2))/(h1*ddz(M)) + &
! & (dhw/dt - dhw0/dt)*0.5*(qwater(M+1,2)-qwater(M,2))
fdiff(1:2) = - lammeth(M)*(qwater(M+1,2,1:2)-qwater(M,2,1:2))/(h1*ddz(M))
! fdiff1 = - diff(ns)*(qsoil(ns)-qsoil(ns-1))/ & ! ml
! & (z(ns)-z(ns-1))
! ===============================================================
!* total methane emission ftot
!* ==============================================================
! if(w.ge.rns)then
! ftot = fdiff + febul + fplant
fdiff_lake_surf(:) = - Flux_atm(:) !- lammeth(1)/(h1*ddz(1))*(qwater(2,2)-qwater(1,2))
! q_sum_new = 0.
! do i = 2, ns-1 !ms+1, ml-1
! q_sum_new = q_sum_new + qsoil(i)*0.5*(z(i+1)-z(i-1))
! end do
! q_sum_new = q_sum_new + qsoil(ns)*0.5*(z(ns)-z(ns-1)) ! ml
! q_sum_new = q_sum_new + qsoil(1)*(0.5*(z(2)-z(1)) + 0.5*ddz(M)*h1)
! do i = 2, ns-1 !ms+1, ml-1
! rprod_sum = rprod_sum + rprod(i)*dt*0.5*(z(i+1)-z(i-1))
!! xx = q_old(i)*0.5*(z(i+1)-z(i-1))+ rprod(i)*dt*0.5*(z(i+1)-z(i-1))
!! oxid_sum = oxid_sum + &
!! & min(oxid(i)*qsoil(i)*0.5*(z(i+1)-z(i-1)), xx)
!! bull_sum = bull_sum + &
!! & min(bull(i)*(qsoil(i)-ch4_crit(i))*dt*0.5*(z(i+1)-z(i-1)), xx) ! cthresh -> ch4_crit(i)
!! plant_sum = plant_sum + &
!! & min(plant(i)*qsoil(i)*dt*0.5*(z(i+1)-z(i-1)), xx)
! oxid_sum = oxid_sum + oxid(i)*qsoil(i)*dt*0.5*(z(i+1)-z(i-1))
! bull_sum = bull_sum + bull(i)*(qsoil(i)-ch4_crit(i))*dt*0.5*(z(i+1)-z(i-1))
! plant_sum = plant_sum + plant(i)*qsoil(i)*dt*0.5*(z(i+1)-z(i-1))
! end do
! rprod_sum = rprod_sum + rprod(ns)*dt*0.5*(z(ns)-z(ns-1)) ! ml
!! xx = q_old(ns)*0.5*(z(ns)-z(ns-1))+ rprod(ns)*dt*0.5*(z(ns)-z(ns-1))
!! oxid_sum = oxid_sum + &
!! & min(oxid(ns)*qsoil(ns)*dt*0.5*(z(ns)-z(ns-1)), xx)
!! bull_sum = bull_sum + &
!! & min(bull(ns)*(qsoil(ns)-ch4_crit(ns))*dt*0.5*(z(ns)-z(ns-1)), xx) ! cthresh -> ch4_crit(ns)
!! plant_sum = plant_sum + &
!! & min(plant(ns)*qsoil(ns)*dt*0.5*(z(ns)-z(ns-1)), xx)
! oxid_sum = oxid_sum + oxid(ns)*qsoil(ns)*dt*0.5*(z(ns)-z(ns-1))
! bull_sum = bull_sum + bull(ns)*(qsoil(ns)-ch4_crit(ns))*dt*0.5*(z(ns)-z(ns-1))
! plant_sum = plant_sum + plant(ns)*qsoil(ns)*dt*0.5*(z(ns)-z(ns-1))
!
! rprod_sum = rprod_sum + rprod(1)*dt*0.5*(z(2)-z(1))
!! xx = q_old(1)*0.5*(z(2)-z(1))+ rprod(1)*dt*0.5*(z(2)-z(1))
!! oxid_sum = oxid_sum + &
!! & min(oxid(1)*qsoil(1)*dt*0.5*(z(2)-z(1)), xx)
!! bull_sum = bull_sum + &
!! & min(bull(1)*(qsoil(1)-ch4_crit(1))*dt*0.5*(z(2)-z(1)), xx) ! cthresh -> ch4_crit(1)
!! plant_sum = plant_sum + &
!! & min(plant(1)*qsoil(1)*dt*0.5*(z(2)-z(1)), xx)
! oxid_sum = oxid_sum + oxid(1)*qsoil(1)*dt*0.5*(z(2)-z(1))
! bull_sum = bull_sum + bull(1)*(qsoil(1)-ch4_crit(1))*dt*0.5*(z(2)-z(1))
! plant_sum = plant_sum + plant(1)*qsoil(1)*dt*0.5*(z(2)-z(1))
! 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.) then
!! Check the conservation of the scheme
! balance = (rprod_sum - oxid_sum - plant_sum - &
! & bull_sum + fdiff*dt - & !fdiff1*dt-
! & (q_sum_new-q_sum_old))/(q_sum_new-q_sum_old)*100. ! procents
! endif
! print*, 'meth_test', qsoil-q_old
! write(*,*) qwater(:,1,1)
! write(*,*) qwater(:,1,2)
deallocate (z)
deallocate (roots)
deallocate (a,b,c,f,y)
return
END SUBROUTINE METHANE2
SUBROUTINE METHANE_OXIDATION2 &
& (M, dt, oxyg, qwater1, qwater2, DIC, Tw)
use METH_OXYG_CONSTANTS
use PHYS_FUNC, only : &
& REACPOT_ARRHEN
implicit none
! Input/output variables
integer(kind=iintegers), intent(in) :: M
real(kind=ireals), intent(in) :: dt
real(kind=ireals), intent(inout) :: oxyg(1:M+1)
real(kind=ireals), intent(inout) :: qwater1(1:M+1)
real(kind=ireals), intent(inout) :: qwater2(1:M+1)
real(kind=ireals), intent(inout) :: DIC(1:M+1)
real(kind=ireals), intent(in) :: Tw(1:M+1)
! Local variables and constants
real(kind=ireals), parameter :: Kelvin0 = 273.15
real(kind=ireals), parameter :: alpha = 0.1
integer(kind=iintegers) :: i, j
real(kind=ireals) :: x(1:2), xx, temp_K ! Work variables
real(kind=ireals) :: qwater_old(1:2), oxyg_old
! write(*,*) qwater1, qwater2
do i = 2, M
temp_K = Tw(i) + Kelvin0
x(1) = dt*REACPOT_ARRHEN(delta_Eq,temp_K,temp0,vq_max) &
& /(k_o2 + oxyg(i))/(k_ch4 + qwater1(i))
x(2) = dt*REACPOT_ARRHEN(delta_Eq,temp_K,temp0,vq_max) &
& /(k_o2 + oxyg(i))/(k_ch4 + qwater2(i))
qwater_old(1) = qwater1(i)
qwater_old(2) = qwater2(i)
oxyg_old = oxyg(i)
do j = 1, 100
! Simple iterations to solve the system of three non-linear equations
xx = oxyg(i)
oxyg(i) = oxyg(i) - &
& alpha*(oxyg(i) - (oxyg_old - 2.*oxyg(i)*(x(1)*qwater1(i) + x(2)*qwater2(i) ) ) )
! qwater1(i) = qwater_old(1) - x(1)*xx*qwater1(i)
qwater1(i) = qwater1(i) - &
& alpha*(qwater1(i) - (qwater_old(1) - x(1)*xx*qwater1(i)))
! qwater2(i) = qwater_old(2) - x(2)*xx*qwater2(i)
qwater2(i) = qwater2(i) - &
& alpha*(qwater2(i) - (qwater_old(2) - x(2)*xx*qwater2(i)))
! write(*,*) oxyg(i) - (oxyg_old - 2.*oxyg(i)*(x(1)*qwater1(i) + x(2)*qwater2(i))), &
! & qwater1(i) - (qwater_old(1) - x(1)*xx*qwater1(i)), &
! & qwater2(i) - (qwater_old(2) - x(2)*xx*qwater2(i)) !, &
! & oxyg(i), qwater1(i), qwater2(i)
enddo
DIC(i) = DIC(i) - 0.5d0*(oxyg(i) - oxyg_old)
enddo
END SUBROUTINE METHANE_OXIDATION2
...@@ -22,7 +22,8 @@ use DRIVING_PARAMS , only : & ...@@ -22,7 +22,8 @@ use DRIVING_PARAMS , only : &
& missing_value, & & missing_value, &
& thermokarst_meth_prod, & ! Switch for old organics methane production for thermokarst lakes & thermokarst_meth_prod, & ! Switch for old organics methane production for thermokarst lakes
& tricemethhydr, & & tricemethhydr, &
& soil_meth_prod ! Switch for methane production from new organics & soil_meth_prod, & ! Switch for methane production from new organics
& soilswitch !switch for simulation of soil processes
use PHYS_CONSTANTS, only : & use PHYS_CONSTANTS, only : &
& row0, g, row0, roa0, roi, & & row0, g, row0, roa0, roi, &
...@@ -125,10 +126,8 @@ implicit none ...@@ -125,10 +126,8 @@ implicit none
real(kind=ireals), intent(in) :: wa(1:gs%ns) ! air content in a soil, kg/kg real(kind=ireals), intent(in) :: wa(1:gs%ns) ! air content in a soil, kg/kg
real(kind=ireals), intent(in) :: Sals(1:gs%ns) ! Salinity, kg/kg real(kind=ireals), intent(in) :: Sals(1:gs%ns) ! Salinity, kg/kg
! COMMON /VEINIT/
real(kind=ireals), intent(in) :: rootss(1:gs%ns) ! meters real(kind=ireals), intent(in) :: rootss(1:gs%ns) ! meters
! COMMON /SOILDAT/
real(kind=ireals), intent(in) :: rhogr(1:gs%ns) ! density of dry soil, kg/m**3 real(kind=ireals), intent(in) :: rhogr(1:gs%ns) ! density of dry soil, kg/m**3
real(kind=ireals), intent(in) :: por(1:gs%ns) ! soil porosity, n/d real(kind=ireals), intent(in) :: por(1:gs%ns) ! soil porosity, n/d
real(kind=ireals), intent(in) :: methgenmh(1:gs%ns) ! The rate of methane generation/consumption due to real(kind=ireals), intent(in) :: methgenmh(1:gs%ns) ! The rate of methane generation/consumption due to
...@@ -190,17 +189,14 @@ implicit none ...@@ -190,17 +189,14 @@ implicit none
real(kind=ireals), save :: fcoarse = 0.5 !the relative volume of the coarse pores real(kind=ireals), save :: fcoarse = 0.5 !the relative volume of the coarse pores
real(kind=ireals), save :: scf = 16.04/(1./864.)! the scale factor mgCH4/m**2/day real(kind=ireals), save :: scf = 16.04/(1./864.)! the scale factor mgCH4/m**2/day
! scf = 12./(24.*1./864)! the scale factor mgC/m**2/h
real(kind=ireals) :: diff(gs%ns),forg(gs%ns),qplant(gs%ns),qebul(gs%ns),q_old(gs%ns),bull(gs%ns),ch4_crit(gs%ns) real(kind=ireals) :: diff(gs%ns),forg(gs%ns),qplant(gs%ns),qebul(gs%ns),q_old(gs%ns),bull(gs%ns),ch4_crit(gs%ns)
real(kind=ireals) :: oxid(gs%ns),dz_soil(gs%ns),z2(gs%ns),fcch4(gs%ns) real(kind=ireals) :: oxid(gs%ns),dz_soil(gs%ns),z2(gs%ns),fcch4(gs%ns)
real(kind=ireals) :: bp(gs%ns), ap(gs%ns) real(kind=ireals) :: bp(gs%ns), ap(gs%ns)
real(kind=ireals) :: rprod(gs%ns), rprod_new(gs%ns), rprod_old(gs%ns), plant(gs%ns) real(kind=ireals) :: rprod(gs%ns), rprod_new(gs%ns), rprod_old(gs%ns), plant(gs%ns)
real(kind=ireals) :: n2(gs%ns) real(kind=ireals) :: n2(gs%ns)
! real(kind=ireals) :: water_porvolrat(1:ns), air_porvolrat(1:ns)
real(kind=ireals), allocatable :: a(:), b(:), c(:), f(:), y(:) real(kind=ireals), allocatable :: a(:), b(:), c(:), f(:), y(:)
! real(kind=ireals), allocatable, save :: mean_rprod_oldc0(:)
real(kind=ireals) :: rnroot real(kind=ireals) :: rnroot
real(kind=ireals) :: froot(gs%ns), anox_crit(gs%ns) real(kind=ireals) :: froot(gs%ns), anox_crit(gs%ns)
...@@ -246,7 +242,12 @@ implicit none ...@@ -246,7 +242,12 @@ implicit none
& f(1:vector_length),y(1:vector_length)) & f(1:vector_length),y(1:vector_length))
a(:) = 0.; b(:) = 0.; c(:) = 0.; f(:) = 0.; y(:) = 0. a(:) = 0.; b(:) = 0.; c(:) = 0.; f(:) = 0.; y(:) = 0.
! dhw = 0.; dhw0 = 0. if (soilswitch%par == 0 .and. botflux_ch4 == missing_value) then
print*, 'Bottom flux of methane (botflux_ch4) must be set if soil is turned off: STOP'
stop
endif
soilswitch_if : if (soilswitch%par == 1) then !Switch for soil processes
z(:) = zsoil(:) z(:) = zsoil(:)
roots(:) = rootss(:) roots(:) = rootss(:)
...@@ -257,7 +258,7 @@ implicit none ...@@ -257,7 +258,7 @@ implicit none
enddo enddo
! TgrAnn(ns) = Tsoil(ns) !ms ! TgrAnn(ns) = Tsoil(ns) !ms
TgrAnn(gs%ns) = 0. ! exludes the factor of annual mean temperature TgrAnn(gs%ns) = 0. ! excludes the factor of annual mean temperature
if (tricemethhydr%par > 0.) then if (tricemethhydr%par > 0.) then
densi => densmh densi => densmh
...@@ -337,17 +338,15 @@ NO3top = 0. ...@@ -337,17 +338,15 @@ NO3top = 0.
xx = HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, Henry_temp_ref, Tsoil(i) + Kelvin0) xx = HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, Henry_temp_ref, Tsoil(i) + Kelvin0)
yy = HENRY_CONST(Henry_const0_n2, Henry_temp_dep_n2, Henry_temp_ref, Tsoil(i) + Kelvin0) yy = HENRY_CONST(Henry_const0_n2, Henry_temp_dep_n2, Henry_temp_ref, Tsoil(i) + Kelvin0)
ch4_crit(i) = rel_conc_ebul_crit*por(i)*(pressoil(i)*xx - xx/yy*n2(i)) ! assuming hydrostatic pressure in soil h1 + z(i) ch4_crit(i) = rel_conc_ebul_crit*por(i)*(pressoil(i)*xx - xx/yy*n2(i)) ! assuming hydrostatic pressure in soil h1 + z(i)
! pressure -> 1.d+5
enddo enddo
do i = 1, gs%ns !ms+1, ml do i = 1, gs%ns
! anox_crit(i) = (por(i)*rhow/rhogr(i)-wi(i)*rhow/rhoi)*0.9
anox_crit(i) = WL_MAX(por(i),rhogr(i),wi(i),tricemethhydr%par)*0.9 anox_crit(i) = WL_MAX(por(i),rhogr(i),wi(i),tricemethhydr%par)*0.9
end do end do
!Calculation of rooting depth !Calculation of rooting depth
rnroot = 0. rnroot = 0.
do i = 1, gs%ns !ms, ml do i = 1, gs%ns
if (roots(i).gt.0.) rnroot = z(i) if (roots(i).gt.0.) rnroot = z(i)
end do end do
...@@ -367,13 +366,6 @@ t_mature = t_grow + 10. ! the temperature at which plants reach maturity, deg. C ...@@ -367,13 +366,6 @@ t_mature = t_grow + 10. ! the temperature at which plants reach maturity, deg. C
!* ============================= !* =============================
anox = 0. anox = 0.
!do i = ns-1, 2, -1 !ml-1, ms+1, -1
! if (z(i) .le. 1.) then ! 100 cm
! if (wl(i) .ge. anox_crit(i)) then
! anox = anox + (z(i+1)-z(i-1))*0.5
! end if
! end if
!end do
if (rnroot > 0.) then if (rnroot > 0.) then
froot(1:gs%ns) = 0. froot(1:gs%ns) = 0.
...@@ -391,19 +383,6 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.) ...@@ -391,19 +383,6 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.)
end do end do
do i = 2, gs%ns do i = 2, gs%ns
! if (wl(i) .lt. anox_crit(i)) then
! diff(i) = diff_unsat*tortuosity_coef*fcoarse ! in the unsaturated
!
! else ! soil layers
! diff(i) = diff_water*tortuosity_coef*fcoarse ! in the water saturated
! diff(i) = DIFF_WATER_METHANE(Tsoil(i)) * &
! & tortuosity_coef*fcoarse
! end if ! soil layers
! diff(i) = 0.25/\cf*( por(i) + por(i-1) )*tortuosity_coef * &
! & ( (water_porvolrat(i) + water_porvolrat(i-1)) * &
! & DIFF_WATER_METHANE( 0.5*(Tsoil(i) + Tsoil(i-1)) ) + &
! & (air_porvolrat(i) + air_porvolrat(i-1)) * &
! & DIFF_AIR_METHANE( 0.5*(Tsoil(i) + Tsoil(i-1)) ) )
diff(i) = tortuosity_coef * DIFF_WATER_METHANE( 0.5*(Tsoil(i) + Tsoil(i-1)) ) diff(i) = tortuosity_coef * DIFF_WATER_METHANE( 0.5*(Tsoil(i) + Tsoil(i-1)) )
end do end do
...@@ -512,7 +491,6 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.) ...@@ -512,7 +491,6 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.)
rprod(i) = rprod_new(i) + rprod_old(i) rprod(i) = rprod_new(i) + rprod_old(i)
! if(qsoil(i).ge.cthresh) then
if (qsoil(i) .ge. ch4_crit(i) .and. wi(i) == 0.e0_ireals) then ! Formation of bubbles is not allowed in frozen soil if (qsoil(i) .ge. ch4_crit(i) .and. wi(i) == 0.e0_ireals) then ! Formation of bubbles is not allowed in frozen soil
fcch4(i) = 1. fcch4(i) = 1.
else else
...@@ -548,7 +526,6 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.) ...@@ -548,7 +526,6 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.)
! methane oxidation occurs only in the unsaturated soil ! methane oxidation occurs only in the unsaturated soil
! layers ! layers
if (wl(i).lt.anox_crit(i)) then if (wl(i).lt.anox_crit(i)) then
! if(z(i).lt.w) then
! oxid(i) = vmax*q10**((Tsoil(i) - TgrAnn(ns))*0.1)/(rkm + qsoil(i)) !ml ! oxid(i) = vmax*q10**((Tsoil(i) - TgrAnn(ns))*0.1)/(rkm + qsoil(i)) !ml
! Currently oxidation is turned off in soil layers. ! Currently oxidation is turned off in soil layers.
...@@ -563,11 +540,7 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.) ...@@ -563,11 +540,7 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.)
! plant(i) = rkp*tveg*froot(i)*fgrow*(1.-0.5)*veg ! plant(i) = rkp*tveg*froot(i)*fgrow*(1.-0.5)*veg
plant(i) = rkp*tveg*froot(i)*fgrow*veg plant(i) = rkp*tveg*froot(i)*fgrow*veg
! plant(i) = 0.
!
! bull(i) - coefficient for ebullition ! bull(i) - coefficient for ebullition
! if(z(i).ge.w) then
if (wl(i) .ge. anox_crit(i)) then if (wl(i) .ge. anox_crit(i)) then
bull(i) = rke*fcch4(i) bull(i) = rke*fcch4(i)
...@@ -601,6 +574,44 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.) ...@@ -601,6 +574,44 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.)
q_old(i) = qsoil(i) q_old(i) = qsoil(i)
end do end do
!=============================================================
! the methane flux due to ebullition (explicit scheme)
!=================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
febul0 = febul0 + 0.5*(z(i+1)-z(i-1))*qebul(i)
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))
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
else
!With soil disactivated, ebullition flux of methane is assumed zero
febul0 = 0.
fbbleflx_ch4_sum = 0.
endif soilswitch_if
! ============================================================= ! =============================================================
! set up the tridiagonal system of equations ! set up the tridiagonal system of equations
! ============================================================= ! =============================================================
...@@ -881,8 +892,6 @@ else ...@@ -881,8 +892,6 @@ else
endif endif
endif endif
! ============================================================ ! ============================================================
! the methane flux due to plant-mediated transport ! the methane flux due to plant-mediated transport
! =================fplant(t)================================== ! =================fplant(t)==================================
...@@ -890,46 +899,23 @@ endif ...@@ -890,46 +899,23 @@ endif
! qplant(i) = rkp*tveg*froot(i)*fgrow* ! qplant(i) = rkp*tveg*froot(i)*fgrow*
! * (1.-0.5)*qsoil(i) ! * (1.-0.5)*qsoil(i)
! end do ! 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. fplant = 0.
if (soilswitch%par == 1) then
do i = 2, gs%ns-1 !ms+1, ml-1 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 fplant = fplant + 0.5*(z(i+1)-z(i-1))*plant(i)*0.5*(qsoil(i) + q_old(i)) * (1 - pox) !* veg
end do end do
endif
ifdeep : if (deepestsoil) then ifdeep : if (deepestsoil) then
! Adding bubble flux from the lowest soil column to horizontally averaged bubble flux !! Adding bubble flux from the lowest soil column to horizontally averaged bubble flux
if (gs%nsoilcols > 1) then !if (gs%nsoilcols > 1) then
call BUBBLEFLUXAVER(gs%ix,gs%iy,gs%nsoilcols) ! call BUBBLEFLUXAVER(gs%ix,gs%iy,gs%nsoilcols)
else !else
fbbleflx_ch4_sum(0) = febul0*fbbleflx_ch4(0) *bathymwater(1) %area_int ! 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(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 ! fbbleflx_ch4_sum(gs%M+1) = febul0*fbbleflx_ch4(gs%M+1)*bathymwater(gs%M+1)%area_int
endif !endif
!The second step of splitting-up scheme - bubble dissolution source in water column !The second step of splitting-up scheme - bubble dissolution source in water column
if (h1 > 0.) then if (h1 > 0.) then
...@@ -992,32 +978,36 @@ endif ifdeep ...@@ -992,32 +978,36 @@ endif ifdeep
! & diff(2)*(qsoil(2) - qsoil(1))*2./(z(2)-z(1)) + & ! & diff(2)*(qsoil(2) - qsoil(1))*2./(z(2)-z(1)) + &
! & xx*( - qebul(1) - plant(1)*qsoil(1) + rprod(1) - oxid(1)*qsoil(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 fdiff1 = - 0.5*diff(gs%ns)*(qsoil(gs%ns) + q_old(gs%ns) - qsoil(gs%ns-1) - q_old(gs%ns-1))/ &
& (z(gs%ns) - z(gs%ns-1)) & (z(gs%ns) - z(gs%ns-1))
else
fdiff1 = 0.
endif
!=============================================================== !===============================================================
! total methane emission ftot ! total methane emission ftot
!============================================================== !==============================================================
!if(w.ge.rns)then
if (deepestsoil) then if (deepestsoil) then
fdiff_lake_surf = Flux_atm !- lammeth(1)/(h1*ddz(1))*(qwater(2,2)-qwater(1,2)) fdiff_lake_surf = Flux_atm !- lammeth(1)/(h1*ddz(1))*(qwater(2,2)-qwater(1,2))
ftot = fdiff_lake_surf + fbbleflx_ch4_sum(0)/bathymwater(1)%area_int + fplant ftot = fdiff_lake_surf + fbbleflx_ch4_sum(0)/bathymwater(1)%area_int + fplant
endif 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
!Diagnostic calculations of methane balance in sediments
q_sum_new = 0. q_sum_new = 0.
do i = 2, gs%ns-1 !ms+1, ml-1 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)) q_sum_new = q_sum_new + qsoil(i)*0.5*(z(i+1)-z(i-1))
end do end do
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(gs%ns)*0.5*(z(gs%ns)-z(gs%ns-1))
q_sum_new = q_sum_new + qsoil(1)*(0.5*(z(2)-z(1))) ! + 0.5*ddz(M)*h1/por(1)) 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 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)) 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)) oxid_sum = oxid_sum + oxid(i)*(qsoil(i) + q_old(i))*dt*0.25*(z(i+1)-z(i-1))
...@@ -1035,7 +1025,7 @@ endif ifdeep ...@@ -1035,7 +1025,7 @@ endif ifdeep
bull_sum = bull_sum + bull(1)*(q_old(1)-ch4_crit(1))*dt*0.5*(z(2)-z(1)) 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)) plant_sum = plant_sum + plant(1)*(qsoil(1) + q_old(1))*dt*0.25*(z(2)-z(1))
! print*, 'bounds_talik', 0.5*(z(i_talik-1)+z(i_talik)), 0.5*(z(i_talik+1)+z(i_talik)) endif
!if( (q_sum_new - q_sum_old).ne.0. .and. gs%isoilcol == gs%nsoilcols) then !if( (q_sum_new - q_sum_old).ne.0. .and. gs%isoilcol == gs%nsoilcols) then
! xx = 0. ! xx = 0.
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment