1.Diagnostics of methane conservation elaborated 2.Misc

source/model/init_var_mod.f90

@@ -251,7 +251,7 @@
   enddo
 
   ! Identical initialization of all soil columns
-  if (nsoilcols > 1) then
+  if (nsoilcols > 1 .and. depth_area(1,2) >= 0.) then
     do i = 1, ns
       wl1   (i,1:nsoilcols-1) = wl1   (i,nsoilcols)
       wi1   (i,1:nsoilcols-1) = wi1   (i,nsoilcols)  

source/model/lake.f90

@@ -2443,6 +2443,7 @@ if (flag_print) then ! The output in ASCII files
     &  work1, z_full, M+1, &
     &  zgrid_out, ngrid_out%par, &
     &  outpath, 'methane_water', i, ndec, .false.)
+    if (nsoilcols > 1) then
       ndec = -1
       i = i + 2
       ! Methane production in soil columns
@@ -2454,6 +2455,7 @@ if (flag_print) then ! The output in ASCII files
       &  rprod_total_newC, work1, nsoilcols, &
       &  work1, nsoilcols, &
       &  outpath, 'rprod_soil', i, ndec, .false.)
+    endif 
     ndec = -1
     i = i + 2
     work1(1:M+1) = carbdi(1:M+1,1)*molm3tomgl_co2 !molm3tonM

source/model/lake_modules.f90

@@ -182,7 +182,7 @@ salice0 = 0.! 4.d-3 ! salinity of ice
 snmeltwat = 0. ! The fraction of snowmelt transferred to water through cracks in ice, 
                ! the rest is frozen on a top of the ice layer
 
-z0_bot = 1.d-3 ! Roughness parameter of bottom surface, m
+z0_bot = 1.d-5 ! Roughness parameter of bottom surface, m
 c_bot = 2.5d-3 ! Momentum exchange coefficient at the sloping bottom
 
 ci_m_roi  = ci * roi

source/model/methane_mod.f90

@@ -589,6 +589,11 @@ forall (i = 1:gs%ns, z(i) > rnroot) forg(i) = veg * exp((rnroot - z(i))*10.)
  q_sum_old = q_sum_old + qsoil(gs%ns)*0.5*(z(gs%ns)-z(gs%ns-1))
  q_sum_old = q_sum_old + 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_old = q_sum_old + qwater(i,1)*h1*ddz05(i-1)
+ enddo
+
  do i = 1, gs%ns !1, ml
    q_old(i) = qsoil(i)
  end do
@@ -925,7 +930,7 @@ ifdeep : if (deepestsoil) then
       qwater(gs%M+1,2) = max(qwater(gs%M+1,2),small_value) !small_value
       qsoil(1) = por(1)*qwater(gs%M+1,2) !porosity added
     elseif (bottom_bc == 2) then
-      qwater(gs%M+1,2) = qwater(gs%M+1,2) + 2.d0*dt/(h1*ddz(gs%M)) & ! porosity added
+      qwater(gs%M+1,2) = qwater(gs%M+1,2) + 2.d0*dt/(h1*ddz(gs%M)) & 
       & *(fbbleflx_ch4_sum(gs%M+1) - fbbleflx_ch4_sum(gs%M)) / bathymwater(gs%M+1)%area_int
       qwater(gs%M+1,2) = max(qwater(gs%M+1,2),small_value) !small_value
     endif
@@ -982,7 +987,6 @@ endif ifdeep
  if (deepestsoil) then 
    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
-   !print*, qwater(:,2), qsoil
  endif
 
  q_sum_new = 0.
@@ -992,6 +996,11 @@ 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))
@@ -1012,9 +1021,19 @@ endif ifdeep
 ! 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 
+ !  xx = 0.
+ !  do i = 1, gs%M+1
+ !    xx = xx + 0.5*(qwater(i,1) + qwater(i,2))*ddz05(i-1)
+ !  enddo
+ !  yy = dhw*(0.5*(qwater(gs%M+1,2)+qwater(gs%M+1,1)) - xx) - &
+ !  & dhw0*0.5*((qwater(gs%M+1,2)+qwater(gs%M+1,1)) - &
+ !  &           (qwater(1     ,2)+qwater(1     ,1)))
  !  !Check the conservation of the scheme in soil
- !  balance = (rprod_sum - oxid_sum - plant_sum - &
- !  &  bull_sum + fdiff*dt - methox2sod*sodbot*dt - & 
+ !  !balance = (rprod_sum - oxid_sum - plant_sum - &
+ !  !&  bull_sum + fdiff*dt - methox2sod*sodbot*dt - & 
+ !  !&  (q_sum_new-q_sum_old)) !/(q_sum_new-q_sum_old)*100. ! percents
+ !  balance = (rprod_sum + yy - oxid_sum - plant_sum - &
+ !  &  bull_sum - Flux_atm*dt - methox2sod*sodbot*dt - & 
  !  &  (q_sum_new-q_sum_old)) !/(q_sum_new-q_sum_old)*100. ! percents
  !  write(*,*) 'balance = ', balance, &
  !  & 'rprod=', rprod_sum, &
@@ -1023,6 +1042,7 @@ endif ifdeep
  !  & 'ebull=', -bull_sum, &
  !  & 'febul0=', -febul0*dt, &
  !  & 'fdiff=',  fdiff*dt, &
+ !  & 'Flux_atm=',  -Flux_atm*dt, &
  !  & 'storage=', - (q_sum_new - q_sum_old)
  !  read*
  !endif