!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*                   SOIL_C/N - Programs                         *!
!*                                                               *!
!*                    Author: F. Suckow                          *!
!*                                                               *!
!*   contains:                                                   *!
!*   SOIL_CN                                                     *!
!*   F_CNV(Cpool, Npool)                                         *!
!*   RMIN_T(temp)                                                *!
!*   RNIT_T(temp)                                                *!
!*   RMIN_W(water, xpv)                                          *!
!*   RNIT_W(water, xpv)                                          *!
!*   RMIN_P(phv)                                                 *!
!*   RNIT_P(phv)                                                 *!
!*   HUMLAY                                                      *!
!*   DECOMP1(Copm, Nopm, cnv, kopm, ksyn, hdiff)                 *!
!*   DECOMP2(Copm, Nopm, cnv, kopm, ksyn, hdiff)                 *!
!*   MINLAY(jlay)                                                *!
!*   N_LEACH(jlay, NH4l, NO3l)                                   *!
!*   S_RESP(Copm_1, Chum_1)                                      *!
!*                                                               *!
!*                  Copyright (C) 1996-2018                      *!
!*     Potsdam Institute for Climate Impact Reserach (PIK)       *!
!*          Authors and contributors see AUTHOR file             *!
!*  This file is part of 4C and is licensed under BSD-2-Clause   *!
!*                   See LICENSE file or under:                  *!
!*     http://www.https://opensource.org/licenses/BSD-2-Clause   *!
!*                           Contact:                            *!
!*       https://gitlab.pik-potsdam.de/foresee/4C                *!
!*                                                               *!
!*****************************************************************!

SUBROUTINE soil_cn

!   Soil C-N budget

use data_climate
use data_out
use data_simul
use data_soil
use data_soil_cn
use data_stand

implicit none

integer j, hnlay, ntr
real Copm_1, Chum_1   ! previous C-content of soil profile
real Nopm_1, Nhum_1   ! previous N-content of soil profile
real Cbc_1, Nbc_1   ! previous C- and N-content of biochar
real Nmin1, N_min_h
type(Coh_Obj), pointer :: p  ! pointer to cohort list

! save previous state of soil C-content
Copm_1 = SUM(C_opm) + C_opm_stem
Chum_1 = SUM(C_hum)
Nopm_1 = SUM(N_opm) + N_opm_stem
Nhum_1 = SUM(N_hum)
N_min_h= N_min
if (flag_bc .gt. 0) then
    Cbc_1  = SUM(C_bc)
    Nbc_1  = SUM(N_bc)
else
    Cbc_1  = 0.
    Nbc_1  = 0.
endif

call humlay                 ! humus layer

! loop over mineral layers
do j=2,nlay
   call minlay(j)
enddo      ! loop over j (nlay)

! soil respiration
call s_resp(Copm_1, Chum_1, Cbc_1)

! daily values
Nleach   = NH4_in + NO3_in
Nupt_d   = SUM(Nupt)
N_an_tot = SUM(NH4) + SUM(NO3)
Nmin1    = Nopm_1 + Nhum_1 - SUM(N_opm) - SUM(N_hum)
if (flag_bc .gt. 0) then
    Nmin1  = Nmin1 + Nbc_1 - SUM(N_bc)
endif  

! yearly cumul. quantities
Nleach_c = Nleach_c + Nleach
Nupt_c   = Nupt_c + Nupt_d
resps_c  = resps_c + respsoil

p => pt%first
do while (associated(p))
    ns = p%coh%species
    ns = p%coh%species
    ntr = p%coh%ntreea
    svar(ns)%Ndem = svar(ns)%Ndem + ntr * p%coh%Ndemc_d
    svar(ns)%Nupt = svar(ns)%Nupt + ntr * p%coh%Nuptc_d

    p%coh%Nuptc_c = p%coh%Nuptc_c + p%coh%Nuptc_d
    p%coh%Ndemc_c = p%coh%Ndemc_c + p%coh%Ndemc_d
    p%coh%N_pool  = p%coh%N_pool  + p%coh%Nuptc_d

    p => p%next
enddo  ! p (cohorts)  

if (flag_dayout .ge. 2) then 
   if (nlay .gt. 6)	then 
      hnlay = 6

   else
      hnlay = nlay
   endif

   N_min_h = N_min - N_min_h
   write (unit_soicna, '(A)') ''
   write (unit_soicnd, '(A)') ''
endif

1000  FORMAT (2I5, 6F10.3, 6F10.1)
1100  FORMAT (2I5, 12F10.3)
1200  FORMAT (2I5, 4F10.3, 4F10.1, F10.2)

END subroutine soil_cn

!**************************************************************

real FUNCTION f_cnv(Cpool, Npool)

! C/N-ratio of a pool
! implicit none

real Cpool, Npool

  if (Npool .lt. 1e-6) then
     f_cnv = 0.
  else     
     f_cnv  = Cpool / Npool
  endif
  
END function f_cnv

!**************************************************************

real FUNCTION rmin_t(temp, rkind)

! reduction of mineralization depending on soil temperature
use data_simul
implicit none

integer rkind
real temp, toptm, Q10

select case (rkind)

case(1)
    toptm  = 35.
    Q10    = 2.9
    rmin_t = exp(log(Q10) * ((temp-toptm)/10.))   ! Stanford

case(2)
    toptm  = 35.
    Q10    = 2.9
    rmin_t = Q10**((temp-toptm)*0.1)              ! van't Hoff 

case(4)
    rmin_t = 1.

case default
    toptm  = 35.
    Q10    = 2.9
    rmin_t = exp(log(Q10) * ((temp-toptm)/10.))   ! Stanford

end select

END function rmin_t

!**************************************************************

real FUNCTION rnit_t(temp, rkind)

! reduction of nitrification depending on soil temperature

implicit none

integer rkind
real temp, toptn, Q10

select case (rkind)

case(1)   ! Stanford
    toptn  = 30.
    Q10    = 2.8
    rnit_t = exp(log(Q10) * ((temp-toptn)/10.))

case(2)   ! van't Hoff
    toptn  = 30.
    Q10    = 2.8
    rnit_t = Q10**((temp-toptn)*0.1)    ! van't Hoff            

case(3)    ! SWAT-approach; Nitrif. only above 5°C
    if (temp .gt. 5.) then
        rnit_t = 0.041 *(temp-5.) 
    else
        rnit_t = 0.
    endif

case(4)
    rnit_t = 1.

case default
    toptn  = 30.
    Q10    = 2.8
    rnit_t = exp(log(Q10) * ((temp-toptn)/10.))   ! Stanford

end select


END function rnit_t

!**************************************************************

real FUNCTION rmin_w(water, xpv)

! reduction of mineralization depending on soil water content
! xpv - pore volume

    rmin_w = 4.0 * water * (1.0-water/xpv) / xpv
if (rmin_w .lt. 0.) rmin_w = 0.

END function rmin_w

!**************************************************************

real FUNCTION rnit_w(water, xpv, xfk, xwp, rkind)

! reduction of nitrification depending on soil water content
! xpv - pore volume

implicit none

integer rkind
real water, xpv, xfk, xwp, nfk, avwat

select case (rkind)

case(1)   ! Franco
    if (water .lt. 0.9*xpv) then
        rnit_w = 4.0 * water * (1.0-water/xpv) / xpv
    else
	    rnit_w = 1.
    endif
    if (rnit_w .lt. 0.) rnit_w = 0.

case(2)  ! SWAT-Ansatz
    nfk = xfk - xwp
    avwat = water - xwp
    if (avwat .lt. 0.25*nfk) then
	    rnit_w = avwat / 0.25 * nfk
    else
	    rnit_w = 1.
    endif

case default
    if (water .lt. 0.9*xpv) then
        rnit_w = 4.0 * water * (1.0-water/xpv) / xpv
    else
	    rnit_w = 1.
    endif
    if (rnit_w .lt. 0.) rnit_w = 0.

end select

END function rnit_w

!**************************************************************
 

real FUNCTION rmin_p(phv)

! reduction of mineralization depending on pH-value
real, dimension(4)  :: a = (/2.5, 4.0, 5.0, 8.0/),  &
                       b = (/0.5, 0.8, 1.0, 1.0/)

call tab_int(a,b,4,phv,value)
rmin_p = value

END function rmin_p

!**************************************************************
 

real FUNCTION rnit_p(phv)

! reduction of nitrification depending on pH-value
real, dimension(4)  :: a = (/2.5, 4.0, 6.0, 8.0/),  &
                       b = (/0.1, 0.3, 1.0, 1.0/)

call tab_int(a,b,4,phv,value)
rnit_p = value

END function rnit_p

!**************************************************************

SUBROUTINE humlay

!   C-N budget of the humus layer
!   (including litter layer)
use data_climate
use data_depo
use data_inter
use data_out
use data_simul
use data_soil
use data_soil_cn
use help_soil_cn
use data_species

implicit none

integer, parameter:: double_prec = kind(0.0D0)
integer i
real (kind = double_prec):: N_hum_1, NH4_1, NO3_1     ! previous state of C- and N-pools
real (kind = double_prec):: N_hum_2, NH4_2, NO3_2     ! actual state of C- and N-pools
real (kind = double_prec):: hnh4, hno3, bilanz, hnhum, hncopm, nh4diff, nhdiff, hdiff, s_hdiff
real (kind = double_prec):: renit			          ! reduction function of nitrif.
real (kind = double_prec):: redtermc, redtermn		  ! red. terms of C-/ N-pools
real Copm, Nopm, hcnv, hcnv_bc, kopm, redopm, Nminl, Nmin1, redbc
logical ldecomp
real, external :: rmin_t, rmin_w, rnit_t, rnit_w, f_cnv 
type (species_litter)  :: sliti

if (flag_dayout .ge. 2) then
  write (unit_soicnr, '(2I5,3E12.3)') time_cur, iday, rmin_t(temps(1), kmint), rmin_w(wats(1), pv(1)), rmin_phv(1)
endif

! reduction factors of mineralization and nitrification
remin = rmin_t(temps(1), kmint) * rmin_w(wats(1), pv(1)) * rmin_phv(1)
renit = rnit_t(temps(1), knitt) * rnit_w(wats(1), pv(1), field_cap(1), wilt_p(1), knitw) * rnit_phv(1)

! add deposition
if (flag_depo .eq. 2) then
    NH_dep = NH_dep * prec_stand    ! conversion g/l in g/m2
    NO_dep = NO_dep * prec_stand
endif
NH4(1) = NH4(1) + NH_dep  
NO3(1) = NO3(1) + NO_dep

Ndep_cum = Ndep_cum + NO_dep + NH_dep

! store state of previous step
N_hum_1  = N_hum(1)
NH4_1    = NH4(1)
NO3_1    = NO3(1)

khr      = k_hum * remin 
hexph    = exp(-khr) 
knr      = k_nit * renit
if (abs(knr-khr) .le. 1E-6) knr = knr + 1E-6
hexpn    = exp(-knr)

! reduction of C- and N-humus-pool by mineralization,
redtermc = C_hum(1) * hexph		  ! part of equation II
redtermn = N_hum_1 * hexph 		  !        -"-

! NH4-pool
if (NH4_1 .gt. 1E-6) then
   term1 = NH4_1 * hexpn 		  ! part of equ. III
else
   term1 = NH4_1
endif
term3    = N_hum_1 * khr * (hexph-hexpn) / (knr-khr) 

if (cnv_hum(1) .lt. 1e-8) cnv_hum(1) = 20.
cnvh     = 1./cnv_hum(1)
redopm   = 1.
redbc    = 1.
slit_1 = slit
ldecomp = .TRUE.
do while (ldecomp)

  ! Decomposition of dead biomass 
   Copm = 0.
   Nopm = 0.
   C_opm_stem = 0.
   N_opm_stem = 0.

   reptermc = 0.
   reptermn = 0.
   term2    = 0.
   term4    = 0.

   s_hdiff = 0.
  ! Decomposition of dead biomass fractions
  do i=1,nspecies
    sliti  = slit_1(i)
    hdiff  = 0.

    if (sliti%C_opm_fol .gt. 1e-8) then 
       kopm = redopm * spar(i)%k_opm_fol
       if (kopm .ge. 1e-8) then
          sliti%cnv_opm_fol  = f_cnv(sliti%C_opm_fol, sliti%N_opm_fol)
          call decomp1(sliti%C_opm_fol, sliti%N_opm_fol, sliti%cnv_opm_fol, &
	           kopm, spar(i)%k_syn_fol, hdiff)   
          s_hdiff = s_hdiff + hdiff
       endif
    endif
   	        
    if (sliti%C_opm_frt(1) .gt. 1e-8) then 
       kopm = redopm * spar(i)%k_opm_frt
       if (kopm .ge. 1e-8) then
          sliti%cnv_opm_frt  = f_cnv(sliti%C_opm_frt(1), sliti%N_opm_frt(1))
          call decomp1(sliti%C_opm_frt(1), sliti%N_opm_frt(1), sliti%cnv_opm_frt, &
	           kopm, spar(i)%k_syn_frt, hdiff) 
          s_hdiff = s_hdiff + hdiff
       endif
    endif
   	        
    if (sliti%C_opm_tb .gt. 1e-8) then 
       kopm = redopm * spar(i)%k_opm_tb
       if (kopm .ge. 1e-8) then
          sliti%cnv_opm_tb  = f_cnv(sliti%C_opm_tb, sliti%N_opm_tb)
          call decomp1(sliti%C_opm_tb, sliti%N_opm_tb, sliti%cnv_opm_tb, &
  	           kopm, spar(i)%k_syn_tb, hdiff) 
          s_hdiff = s_hdiff + hdiff
       endif
    endif

    select case (flag_decomp)
    case (0, 10, 20, 30, 40)
       if (sliti%C_opm_crt(1) .gt. 1e-8) then 
         kopm = redopm * spar(i)%k_opm_crt
         if (kopm .ge. 1e-8) then
            sliti%cnv_opm_crt  = f_cnv(sliti%C_opm_crt(1), sliti%N_opm_crt(1))
            call decomp1(sliti%C_opm_crt(1), sliti%N_opm_crt(1), sliti%cnv_opm_crt, &
	             kopm, spar(i)%k_syn_crt, hdiff) 
            s_hdiff = s_hdiff + hdiff
         endif
       endif
   	        
       if (sliti%C_opm_stem .gt. 1e-8) then 
         kopm = redopm * spar(i)%k_opm_stem
         if (kopm .ge. 1e-8) then
            sliti%cnv_opm_stem  = f_cnv(sliti%C_opm_stem, sliti%N_opm_stem)
            call decomp1(sliti%C_opm_stem, sliti%N_opm_stem, sliti%cnv_opm_stem, &
	              kopm, spar(i)%k_syn_stem, hdiff) 
            s_hdiff = s_hdiff + hdiff
         endif
       endif

   case (1, 11, 21, 31, 41)
       if (sliti%C_opm_crt(1) .gt. 1e-8) then 
         kopm = redopm * spar(i)%k_opm_crt
         if (kopm .ge. 1e-8) then
            sliti%cnv_opm_crt  = f_cnv(sliti%C_opm_crt(1), sliti%N_opm_crt(1))
            call decomp2(sliti%C_opm_crt(1), sliti%N_opm_crt(1), sliti%cnv_opm_crt, &
	             kopm, spar(i)%k_syn_crt, hdiff) 
            s_hdiff = s_hdiff + hdiff
         endif
       endif
   	        
       if (sliti%C_opm_stem .gt. 1e-8) then 
         kopm = redopm * spar(i)%k_opm_stem
         if (kopm .ge. 1e-8) then
            sliti%cnv_opm_stem  = f_cnv(sliti%C_opm_stem, sliti%N_opm_stem)
            call decomp2(sliti%C_opm_stem, sliti%N_opm_stem, sliti%cnv_opm_stem, &
	              kopm, spar(i)%k_syn_stem, hdiff) 
            s_hdiff = s_hdiff + hdiff
         endif
       endif
  
   end select

  ! pools of dead biomass without stems
   Copm = Copm + sliti%C_opm_fol + sliti%C_opm_frt(1) + sliti%C_opm_crt(1) + sliti%C_opm_tb
   Nopm = Nopm + sliti%N_opm_fol + sliti%N_opm_frt(1) + sliti%N_opm_crt(1) + sliti%N_opm_tb

  ! dead stems
   C_opm_stem = C_opm_stem + sliti%C_opm_stem
   N_opm_stem = N_opm_stem + sliti%N_opm_stem

   slit(i) = sliti   	        
  
  enddo 
  
! Decomposition of biochar
  if (flag_bc .gt. 0) then
    if (C_bc(1) .gt. 1e-8) then 
       kbc = redbc * k_bc
       if (kbc .ge. 1e-8) then
          hcnv_bc  = f_cnv(C_bc(1), N_bc(1))
          call decomp1(C_bc(1), N_bc(1), hcnv_bc, kbc, k_syn_bc, hdiff) 
          s_hdiff = s_hdiff + hdiff
       endif
    endif
  endif

  
  ldecomp = .FALSE.

  C_opm(1) = Copm
  N_opm(1) = Nopm

! C- and N-humus-pool: reduction by mineralization, supply by turnover of organic primary matter
  C_hum(1) = redtermc + reptermc
  N_hum_2  = redtermn + reptermn
  N_hum(1) = N_hum_2

! ammonium pool
  hnh4    = term1 + term2 + term3 + khr/(knr-khr) * term4
  NH4(1)  = hnh4
  nhdiff  = N_hum_1 - N_hum_2 
  nh4diff = NH4_1 - NH4(1)
  Nminl = hnh4 - NH4_1 - NO3(1)    ! daily net min.

! nitrat pool from balance
  hno3   = NO3_1 + s_hdiff + nhdiff + nh4diff
  NO3(1) = hno3

  if (hnh4 .lt. 0.0 .or. hno3 .lt. 0.0)   then
     redopm = 0.9 * redopm
     if (redopm .ge.  1E-8) then
        ldecomp = .TRUE.
     else
        if (NH4(1) .lt. 1E-10) NH4(1) = 0.
        if (NO3(1) .lt. 1E-10) NO3(1) = 0.
     endif
  endif

  Nminl = Nminl + NO3(1)      ! daily net min. per layer
enddo   ! ldecomp

Nmin(1) = Nminl
N_min = N_min + Nminl      ! cumul. yearly net min.

call n_leach(1)	! without balance

! new balance after leaching
NH4(1) = NH4(1) - NH4_in
NO3(1) = NO3(1) - NO3_in

call n_upt(1)	! with balance

if (flag_dayout .ge. 2) then
  write (unit_soicna, '(2I5,E12.3)', advance='no') time_cur, iday, remin
  write (unit_soicnd, '(2I5,E12.3)', advance='no') time_cur, iday, Nminl
endif

END subroutine humlay

!**************************************************************

SUBROUTINE decomp1(Copm, Nopm, cnv, kopm, ksyn, hdiff)

! Decomposition of dead biomass fractions per species

use help_soil_cn

implicit none

integer, parameter:: double_prec = kind(0.0D0)
real Copm, Nopm        ! C- and N-pool of primary organic matter fraction
real kopm, ksyn        ! mineralisation and synthesis coeff. of opm-fraction 
real kor               ! reduced mineralisation coeff. of opm-fraction 
real N_opm_1, C_opm_1  ! previous state of C- and N-pools
real hexpo             ! exponential part
real cnv               ! C/N-ratio of opm-fraction
real exterm
real (kind = double_prec):: hdiff 
real gamma 

   ! store state of previous step
   C_opm_1  = Copm     
   N_opm_1  = Nopm	   

   kor = kopm * remin	   ! reduction of miner. coeff.
  ! avoid denominators near zero
   if (abs(kor-khr) .lt. 1E-6) kor = kor + 1E-6
   if (abs(kor-knr) .lt. 1E-6) kor = kor + 1E-6
   hexpo    = exp(-kor)
   Copm = C_opm_1 * hexpo	   ! equations II
   Nopm = N_opm_1 * hexpo	   !     -"-

   ! reproduction of C- and N-humus-pool by turnover of organic primary matter
   exterm   = hexph - hexpo
   gamma    = cnv * cnvh

  if (abs(kor-khr) .gt. 1E-6) then
   reptermc = reptermc + C_opm_1 * ksyn * kor * exterm / (kor-khr)      ! part of equ. II
   reptermn = reptermn + N_opm_1 * gamma*ksyn * kor * exterm / (kor-khr)  ! part of equ. II
  endif

  ! change of ammonium pool
  if (abs(kor-knr) .gt. 1E-6) then
   term2 = term2 + (1.-gamma*ksyn)*kor * N_opm_1 * (hexpn - hexpo) / (kor - knr)	  ! part of equ. III
  endif
  if ((abs(kor-khr) .gt. 1E-6) .and. (abs(kor-knr) .gt. 1E-6)) then
   term4 = term4 + gamma*ksyn*kor * N_opm_1 									    & ! part of equ. III
                 * ((kor-khr) * hexpn + (knr-kor) * hexph + (khr-knr) * hexpo)  &
		         / ((khr - kor) * (kor - knr))
  endif
   
   hdiff = N_opm_1 - Nopm   ! N-change rate in organic primary matter 
   
END subroutine decomp1

!**************************************************************

SUBROUTINE decomp2(Copm, Nopm, cnv, kopm, ksyn, hdiffn)

! Decomposition of dead stem biomass per species

use help_soil_cn

implicit none

integer, parameter:: double_prec = kind(0.0D0)
real Copm, Nopm        ! C- and N-pool of primary organic matter fraction
real kopm, ksyn        ! mineralisation and synthesis coeff. of opm-fraction 
real kor               ! reduced mineralisation coeff. of opm-fraction 
real N_opm_1, C_opm_1  ! previous state of C- and N-pools
real hexpo             ! exponential part
real cnv               ! C/N-ratio of opm-fraction
real (kind = double_prec):: hdiffn, hdiffc 

   ! store state of previous step
   C_opm_1  = Copm     
   N_opm_1  = Nopm	   

   kor = kopm * remin	   ! reduction of miner. coeff.
!    avoid denominators near zero
   if (abs(kor) .lt. 1E-6) kor = kor + 1E-6
   hexpo    = exp(-kor)
   Copm = C_opm_1 * hexpo	   ! equations II
   Nopm = N_opm_1 * hexpo	   !     -"-

   ! reproduction of C- and N-humus-pool by turnover of organic primary matter
   hdiffn = N_opm_1 - Nopm   ! N-change rate in organic primary matter 
   hdiffc = hdiffn / cnvh
   reptermn = reptermn + hdiffn
   reptermc = reptermc + hdiffc 
  
END subroutine decomp2

!**************************************************************

SUBROUTINE minlay(jlay)

!   C-N budget of a mineral layer

use data_climate
use data_out
use data_simul
use data_soil
use data_soil_cn
use help_soil_cn
use data_species

implicit none

! input:
integer jlay         ! number of actual layer

!------------------------------------------------------------

integer, parameter:: double_prec = kind(0.0D0)
integer i
real (kind = double_prec):: N_hum_1, NH4_1, NO3_1     ! previous state of C- and N-pools
real (kind = double_prec):: hnh4, hno3, bilanz, hnhum, hncopm, nh4diff, nhdiff, hdiff, s_hdiff
real (kind = double_prec):: renit			          ! reduction function of nitrif.
real (kind = double_prec):: redtermc, redtermn		  ! red. terms of C-/ N-pools
real Copm, Nopm, hcnv, hcnv_bc, kopm, redopm, Nminl, Nmin1, redbc
real, dimension(nspecies):: Copm_frt_1, Nopm_frt_1, Copm_crt_1, Nopm_crt_1
logical ldecomp
real, external :: rmin_t, rmin_w, rnit_t, rnit_w, f_cnv 

! reduction factors of mineralization and nitrification
remin = rmin_t(temps(jlay), kmint) * rmin_w(wats(jlay), pv(jlay)) * rmin_phv(jlay)
renit = rnit_t(temps(jlay), knitt)  * rnit_phv(jlay) * &
        rnit_w(wats(jlay), pv(jlay), field_cap(jlay), wilt_p(jlay), knitw)

if (flag_dayout .eq. 3) then
    write (1122, *) 'minlay ', iday, jlay
endif    

! add N transport from above layer
NH4(jlay) = NH4(jlay) + NH4_in
NO3(jlay) = NO3(jlay) + NO3_in

! store state of previous step
N_hum_1  = N_hum(jlay)
NH4_1    = NH4(jlay)
NO3_1    = NO3(jlay)
Nopm_frt_1   = slit%N_opm_frt(jlay)
Copm_frt_1   = slit%C_opm_frt(jlay)
Nopm_crt_1   = slit%N_opm_crt(jlay)
Copm_crt_1   = slit%C_opm_crt(jlay)
redopm   = 1.
redbc    = 1.

khr      = k_hum_r * remin 
hexph    = exp(-khr) 
knr      = k_nit * renit
if (abs(knr-khr) .le. 1E-6) knr = knr + 1E-6
hexpn    = exp(-knr)

! reduction of C- and N-humus-pool by mineralization,
redtermc = C_hum(jlay) * hexph          ! part of equation II
redtermn = N_hum_1 * hexph              !        -"-

! NH4-pool
term1    = NH4_1 * hexpn 		  ! part of equ. III
term3    = N_hum_1 * khr * (hexph-hexpn) / (knr-khr) 

if (cnv_hum(jlay) .lt. 1e-8) then
  if (cnv_hum(jlay-1) .ge. 1e-8) then  
      cnv_hum(jlay) = cnv_hum(jlay-1) 
  else    
      cnv_hum(jlay) = 20.
  endif
endif
cnvh = 1./cnv_hum(jlay) 

ldecomp = .TRUE.
do while (ldecomp)
  ! Decomposition of dead biomass 
   reptermc = 0.
   reptermn = 0.
   term2    = 0.
   term4    = 0.
   s_hdiff = 0.
   do i=1,nspecies
      Nopm = Nopm_frt_1(i)
      kopm = redopm * spar(i)%k_opm_frt
      if (Nopm .ge. 1e-8 .and. kopm .ge. 1e-8) then
         Copm = Copm_frt_1(i)
         hcnv = f_cnv(Copm, Nopm)
          if ((time .eq.1) .and. (jlay .gt. 155)) then
          endif    
         call decomp1(Copm, Nopm, hcnv, kopm, spar(i)%k_syn_frt, hdiff) 
      
         slit(i)%C_opm_frt(jlay) = Copm
         slit(i)%N_opm_frt(jlay) = Nopm
         cnv_opm(jlay) = hcnv
      else
         hdiff = 0.
      endif   ! Nopm
      s_hdiff = s_hdiff + hdiff

      Nopm = Nopm_crt_1(i)
      kopm = redopm * spar(i)%k_opm_crt
      if (Nopm .ge. 1e-8 .and. kopm .ge. 1e-8) then
         Copm = Copm_crt_1(i)
         hcnv = f_cnv(Copm, Nopm)
          if ((time .eq.1) .and. (jlay .gt. 155)) then
          endif    
         select case (flag_decomp)
         case (0, 10, 20, 30, 40)
            call decomp1(Copm, Nopm, hcnv, kopm, spar(i)%k_syn_crt, hdiff) 
      
         case (1, 11, 21, 31, 41)
            call decomp2(Copm, Nopm, hcnv, kopm, spar(i)%k_syn_crt, hdiff) 
         end select
         slit(i)%C_opm_crt(jlay) = Copm
         slit(i)%N_opm_crt(jlay) = Nopm
         cnv_opm(jlay) = hcnv
      else
         hdiff = 0.
      endif   ! Nopm
      s_hdiff = s_hdiff + hdiff

   enddo   ! nspecies
  
   ! Decomposition of biochar
   if (flag_bc .gt. 0) then
     if (C_bc(jlay) .gt. 1e-8) then 
       kbc = redbc * k_bc
       if (kbc .ge. 1e-8) then
          hcnv_bc  = f_cnv(C_bc(jlay), N_bc(jlay))
          call decomp1(C_bc(jlay), N_bc(jlay), hcnv_bc, kbc, k_syn_bc, hdiff) 
          s_hdiff = s_hdiff + hdiff
       endif
     endif
   endif

   ldecomp = .FALSE.

   C_opm(jlay) = SUM(slit%C_opm_frt(jlay)) + SUM(slit%C_opm_crt(jlay))
   N_opm(jlay) = SUM(slit%N_opm_frt(jlay)) + SUM(slit%N_opm_crt(jlay))

  ! C- and N-humus-pool: reduction by mineralization,
  !                      supply by turnover of organic primary matter
   C_hum(jlay) = redtermc + reptermc
   hnhum       = redtermn + reptermn
   N_hum(jlay) = hnhum

   ! ammonium pool
   hnh4      = term1 + term2 + term3 + khr/(knr-khr) * term4
   NH4(jlay) = hnh4
   nhdiff    = N_hum_1 - N_hum(jlay) 
   nh4diff   = NH4_1 - NH4(jlay)
   bilanz    = NO3(jlay) + s_hdiff   &
               + nhdiff + nh4diff
   Nminl = NH4(jlay) - NH4_1 - NO3(jlay)     ! daily net min.

   ! nitrate pool from balance
   hno3 = NO3_1 + s_hdiff + nhdiff + nh4diff
   NO3(jlay) = hno3

   if (hnh4 .lt. 0.0 .or. hno3 .lt. 0.0)   then
       redopm = 0.9 * redopm
      if (redopm .ge.  1E-8) then
         ldecomp = .TRUE.
      else
         if (NH4(jlay) .lt. 1E-10) NH4(jlay) = 0.
         if (NO3(jlay) .lt. 1E-10) NO3(jlay) = 0.
      endif
   endif

      Nminl   = Nminl + NO3(jlay)     ! daily net min. per layer
      bilanz = bilanz - NO3(jlay)
enddo   ! ldecomp

Nmin(jlay) = Nminl
N_min = N_min + Nminl        ! cumul. yearly net min.

call n_leach(jlay)	! without balance

! new balance after leaching
NH4(jlay) = NH4(jlay) - NH4_in
NO3(jlay) = NO3(jlay) - NO3_in

call n_upt(jlay)	! with balance

if (flag_dayout .ge. 2) then
  write (unit_soicna, '(E12.3)', advance='no') remin
  write (unit_soicnd, '(E12.3)', advance='no') Nminl
endif

END subroutine minlay

!**************************************************************

SUBROUTINE n_leach(jlay)

!   N leaching and new balance
!   Addition of deposition to the anorganic pools

use data_climate
use data_simul
use data_soil
use data_soil_cn
use help_soil_cn
use data_species

implicit none

! input:
integer jlay        ! number of actual layer

!-----------------------------------------------------------

real NH4f, NO3f  ! free available NH4-, NO3-N
real perc_w      ! relative part of percolated water

! NH4 and NO3 partly fixed

if (NH4(jlay) .lt. 1E-25) then
continue
endif

NH4f   = NH4(jlay) * pNH4f
NO3f   = NO3(jlay) * pNO3f

! relative part of percolated water
perc_w = perc(jlay) / (wats(jlay) + perc(jlay) + wupt_r(jlay) + wupt_ev(jlay))

! N transport
NH4_in = NH4f * perc_w
NO3_in = NO3f * perc_w

END subroutine n_leach

!**************************************************************


SUBROUTINE s_resp(Copm_1, Chum_1, Cbc_1)

! Estimation of soil respiration   

use data_climate
use data_simul
use data_soil
use data_soil_cn
use help_soil_cn
use data_species

implicit none

! input:
real Copm_1, Chum_1, Cbc_1   ! previous C-content of soil profile
real Sum_C_opm, Sum_C_hum, Sum_C_bc

!-----------------------------------------------------------

Sum_C_opm = SUM(C_opm) + C_opm_stem
Sum_C_hum = SUM(C_hum)
respsoil = Copm_1 + Chum_1 - Sum_C_opm - Sum_C_hum  
if (flag_bc .gt. 0) then
    Sum_C_bc = SUM(C_bc)
    respsoil = respsoil + Cbc_1 - Sum_C_bc  
endif

END subroutine s_resp

!**************************************************************

SUBROUTINE s

!   

use data_climate
use data_simul
use data_soil
use data_soil_cn
use help_soil_cn
use data_species

implicit none


END subroutine s

!**************************************************************




 
 
