!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*                Soil and Water - Programs                      *!
!*                                                               *!
!*   contains:                                                   *!
!*   SOIL                                                        *!
!*   SOIL_INI                                                    *!    
!*   SOIL_WAT                                                    *!
!*   UPT_WAT                                                     *!
!*   FRED1 - ...11                                                *!
!*   TAKE_WAT                                                    *!
!*   BUCKET                                                      *!
!*   SNOWPACK                                                    *!
!*   HUM_ADD                                                     *!
!*   BC_APPL: application of biochar                             *!
!*                                                               *!
!*                  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

!   Soil processes (frame)

use data_climate
use data_depo
use data_out 
use data_simul
use data_soil
use data_soil_cn

implicit none

call evapo
call intercep
call soil_wat
call soil_temp

if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
    call soil_stress	!calculate ground stress factors
    call cr_depth		!define root depth 
endif
call soil_cn
call root_eff

END subroutine soil

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

SUBROUTINE soil_ini

! Initialisation of soil data and parameters

use data_inter
use data_evapo
use data_out
use data_par
use data_simul
use data_soil
use data_soil_cn
use data_species
use data_stand

implicit none

integer i,j,k
real d_0, t_0
! Table of quarz and clay content (mass%) versus wlam
real, dimension(17) :: xwlam  = (/1.5, 1.15, 0.9, 0.67, 0.6,  0.5, 0.38, 0.37, 0.3, 0.29, 0.27, 0.26, 0.25, 0.24, 0.23, 0.22, 0.15/),  &
                       yquarz = (/93.0,85.0,80.0, 82.0,76.0, 64.0, 65.0, 51.0,60.0, 30.0, 14.0, 10.0, 12.0, 20.0, 30.0, 43.0, 23.0/),  &
                       yclay  = (/3.0, 3.0,  3.0, 12.0, 6.0,  6.0, 10.0,  4.0,21.0, 12.0, 10.0, 37.0, 15.0, 40.0, 30.0, 35.0, 55.0/) 
real value
real, dimension(nlay):: humush(nlay)
real, dimension(nlay):: xfcap, xwiltp, xpv   ! output of addition for water capacities

! estimation of soil layer values
d_0 = 0.
do j = 1, nlay
   t_0      = thick(j)
   mid(j)   = d_0 + 0.5*t_0
   d_0      = d_0 + t_0
   depth(j) = d_0
enddo

perc	= 0.
wupt_r	= 0.
wupt_ev	= 0.
thick_1 = thick(1)

select case (flag_soilin)
case (0,2)
    do i=1,nlay
        if (i .gt. 1) then
            call tab_int(xwlam, yquarz,  17, wlam(i), value)
            sandv(i)  = value / 100.   ! Mass% of mineral fraction
            call tab_int(xwlam, yclay,  17, wlam(i), value)
            clayv(i)  = value / 100.
            siltv(i)  = 1. - clayv(i) - sandv(i)
        else
            sandv(1)  = 0.0
            clayv(1)  = 0.0
            siltv(1)  = 0.0
        endif
    enddo

case (1,3,4)
    clayv  = clayv / 100.
    sandv  = sandv / 100.
    siltv  = 1. - clayv - sandv
    if ((sandv(1) .le. zero) .and. (clayv(1) .le. zero)) siltv(1) = 0.
    skelv  = skelv / 100.
    humusv = humusv / 100.
end select

! Settings for subroutine take_wat
skelfact  = 1.
pv        = skelfact * pv_v * thick * 0.1	       ! mm
wilt_p    = skelfact * wilt_p_v * thick * 0.1	   ! mm
field_cap = skelfact * f_cap_v * thick * 0.1	   ! mm
wats      = field_cap   ! mm
watvol    = f_cap_v                    ! vol%

n_ev_d = nlay
nlgrw  = nlay+1
do i=1,nlay
	if (w_ev_d .gt. depth(i)) n_ev_d = i
	if (grwlev .gt. depth(i)) nlgrw = i+1
    vol(i)    = thick(i) * 10000.  
enddo

! dry mass of first layer
dmass = vol * dens
rmass1 = dmass(1) - (C_hum(1) + C_opm(1)) / cpart    ! corection term of first layer

humush = humusv
if (2.*C_hum(1) .lt. humusv(1)*dmass(1)) then
    humusv(1) = C_hum(1) / (dmass(1) * cpart)
endif
do i=2, nlay
    humusv(i) = C_hum(i) / (dmass(i) * cpart)
enddo

if (flag_bc .gt. 0) y_bc_n = 1      ! actual number of biochar application

! calculation of additions for water capacities 
call hum_add(xfcap, xwiltp, xpv)

fcaph  = f_cap_v - xfcap
wiltph = wilt_p_v - xwiltp
pvh    = pv_v - xpv

! ground water
do i = nlgrw, nlay
    wats(i) = pv(i)
enddo

interc_can  = 0.
int_st_can  = 0.
interc_sveg = 0.
int_st_sveg = 0.
aev_i       = 0.




wat_tot = SUM(wats)
END subroutine soil_ini

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

SUBROUTINE soil_wat
use data_out
! soil water balance

use data_climate
use data_evapo
use data_inter
use data_out
use data_par 
use data_simul
use data_soil
use data_species
use data_stand

implicit none
 
real :: eva_dem		! evaporation demand of soil
real :: p_inf = 0.  ! infiltrated water
real :: pev		    ! local: soil evaporation
real :: watot, wetot    ! total water content at start and end
real :: wutot           ! total water uptake from the soil
real :: wutot_ev        ! total water uptake by soil evaporation
real :: wutot_r         ! total water uptake by roots
real, allocatable, dimension(:) :: upt   ! local array for uptake
real, external :: wat_new

real enr, wa, we, percolnl, snow_sm, buckdepth
integer j

allocate (upt(nlay))
wupt_ev = 0.
aev_s = 0.
prec_stand = MAX(0., prec - interc_can - interc_sveg)    ! stand precipitation	

if (flag_int .gt. 1000) then
    prec_stand = prec_stand * prec_stand_red / 100.
endif

call snowpack(snow_sm, p_inf, pev)

if (anz_coh .le. 0) pev = pet
eva_dem = MAX(0., p_inf) - pev		! evaporation demand of soil

!     if all stand precipitation is evaporated and there is still a demand
!     there is an uptake from soil layers (up to an certain depth)
if (eva_dem .lt. 0.) then
   if (snow .le. 0) call take_wat(eva_dem, cover) 
   aev_s = aev_s + p_inf + SUM(wupt_ev)
else
   aev_s = aev_s + pev
endif	! eva_dem

aet   = aev_s + aev_i
p_inf = MAX(eva_dem, 0.)
upt   = wupt_ev
watot = SUM(wats)			! total initial water content

do j = 1, nlgrw-1

  enr   = p_inf - upt(j)
  wa    = wats(j) - field_cap(j)
  we    = wat_new(wa, enr, j)
  p_inf = enr + wa - we

  perc(j) = MAX(p_inf, 0.)
  wats(j) = MAX(we+field_cap(j), wilt_p(j))

enddo

do j = nlgrw, nlay

  enr   = p_inf - upt(j)
  wa    = wats(j) - field_cap(j)
  we    = pv(j) - field_cap(j)   ! ground water level is constant!
  p_inf = enr + wa - we
  perc(j) = MAX(p_inf, 0.)
  wats(j) = MAX(we+field_cap(j), wilt_p(j))

enddo

if (flag_wred .le. 10) then
   call upt_wat
else
   call upt_wat1
endif

! root uptake balanced imediate after calculation
upt  = upt + wupt_r
wats = wats - wupt_r

watvol = 10.*wats/(thick * skelfact)   ! estimation for complete layer in Vol% without skeleton (only soil substrate)

! total water quantities
wetot    = SUM(wats)		! total final water content
wutot_ev = SUM(wupt_ev)		! total water uptake by soil evaporation
wutot_r  = SUM(wupt_r)		! total water uptake by roots
wutot    = wutot_ev + wutot_r	! total water uptake
aet      = aet + wutot_r        ! daily total aet

percolnl = perc(nlay)

trans_tree = 0.
trans_sveg = 0.
   zeig => pt%first
   do while (associated(zeig))
     if (zeig%coh%species .le. nspec_tree) then
        trans_tree = trans_tree + zeig%coh%supply 
     else
        trans_sveg = trans_sveg + zeig%coh%supply 
     endif
     zeig => zeig%next
   enddo  ! zeig (cohorts) 

! cumulative water quantities
perc_cum   = perc_cum + perc(nlay)
wupt_r_c   = wupt_r_c + wutot_r
wupt_e_c   = wupt_e_c + wutot_ev
wupt_cum   = wupt_cum + wutot
aet_cum    = aet_cum + aet
dew_cum    = dew_cum + dew_rime
tra_tr_cum = tra_tr_cum + trans_tree 
tra_sv_cum = tra_sv_cum + trans_sveg

call bucket(bucks_100, bucks_root, buckdepth)

! number of drought days per layer
where ((wats-0.2) .le. wilt_p) s_drought = s_drought+1

 if (flag_dayout .ge. 2) then
   write (unit_wat, '(2I5, 7F7.2, 24F8.2)') time_cur, iday, airtemp, prec, interc_can, int_st_can, &
                    interc_sveg, int_st_sveg, snow, snow_sm,  pet, trans_dem, &
                    pev, aev_s, aev_i, perc(nlay), watot, wetot, wutot, wutot_ev, wutot_r,&
                    trans_tree,trans_sveg, eva_dem, gp_can, aet, ceppot_can, ceppot_sveg
endif
deallocate (upt)
		    
END	subroutine soil_wat

!**************************************************************
SUBROUTINE upt_wat

! Water uptake by roots

use data_simul
use data_evapo
use data_soil
use data_stand 
use data_par
use data_species
use data_climate

implicit none

real, dimension(1:anz_coh) :: tr_dem   ! auxiliary arrays for cohorts
real wat_ava, hdem, frtrel, frtrel_1, hupt, hupt_c, totfrt_2, hv, demand_mistletoe_canopy
real wat_at    ! total available water per layer
real wat_ar    ! total available water per layer with uptake resistance
real hred      ! resistance coefficient
real, external :: fred1, fred2, fred3, fred4, fred5, fred6, fred7, fred11
integer i, ianz, j, nroot3


! Calculation of Water Demand
ianz = anz_coh

tr_dem = 0
hdem   = 0
hv     = pet-aev_i-pev_s
if (hv .lt. 0.) hv = 0.

select case (flag_eva)
case (0,1,3)
 if((pet .gt. 0.) .and. (hv .gt. 0.)) then
    trans_dem = hv * alfm * (1. - exp(-gp_tot/gpmax))  ! pet (potential evapotranspiration) is reduced by intereption evaporation and potential ground evaporation
 else
 trans_dem = 0.0 
 endif 
   if (gp_tot .gt. zero) then
      hdem   = trans_dem / gp_tot
   else
      hdem= 0.
   endif
case (8)  
 ! potential transpiration demand of each cohort
   if ((gp_tot .gt. zero) .and. (hv .gt. 0.)) then
      hdem   = (pet-aev_i-aev_s) / gp_tot
   else
      hdem= 0.
   endif

case (2,4)
   trans_dem = 0.

case (6,16,36)
 ! Eucalyptus
 hv = pet

 if((pet .gt. 0.) .and. (hv .gt. 0.)) then
    trans_dem = hv * alfm * (1. - exp(-gp_tot/gpmax)) 
 else
    trans_dem = 0.
 endif 

 ! preparation: potential transpiration demand of each cohort
   if (gp_tot .gt. zero) then
      hdem   = trans_dem / gp_tot
   else
      hdem= 0.
   endif

case (7,17,37)
 trans_dem = hv

 ! potential transpiration demand of each cohort
   if (gp_tot .gt. zero) then
      hdem   = trans_dem / gp_tot
   else
      hdem= 0.
   endif

end select

hdem = max(0., hdem)

! Distribution of total Demand into Demands of Cohorts

!extraction of demand of mistletoe cohort (for case flag eva = 1,3,6,7...)
zeig => pt%first
 do while (associated(zeig))
  if (zeig%coh%species.eq.nspec_tree+2) then
      demand_mistletoe_canopy=zeig%coh%gp * zeig%coh%ntreea * hdem
  end if
  zeig => zeig%next
 enddo
zeig => pt%first
i = 1
do while (associated(zeig))

  select case (flag_eva)
  case (0, 1, 3, 6, 7, 16, 17, 36, 37)

            !uppermost tree cohort (with flag mistletoe) gets additinal demand of mistletoe
            if (zeig%coh%mistletoe.eq.1) then
                zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem + demand_mistletoe_canopy
            elseif (zeig%coh%species.eq.nspec_tree+2) then                ! set demand of mistletoe to zero as it will be fullfilled by the tree
               zeig%coh%demand=0.                                         ! set to zero because demand has been added to the infested tree cohort
            else
                zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem    ! all other cohorts get their demand
            end if

  case (2,4)
 !uppermost tree cohort (with flag mistletoe) gets additinal demand, that of mistletoe
            if (zeig%coh%mistletoe.eq.1) then
                zeig%coh%demand = (max(0., zeig%coh%demand - zeig%coh%aev_i)  + demand_mistletoe_cohort)
            endif
            if (zeig%coh%species.eq.nspec_tree+2) then                ! set demand of mistletoe to zero as it will be fullfilled by the tree
                zeig%coh%demand=0.
            endif
            if (zeig%coh%mistletoe.ne.1 .AND. zeig%coh%species.ne.nspec_tree+2) then
               zeig%coh%demand = max(0., zeig%coh%demand - zeig%coh%aev_i)   
            end if
            trans_dem = trans_dem + zeig%coh%demand
    end select

  tr_dem(i) = zeig%coh%demand   ! demand of transpiration per cohort 
   i = i + 1
  zeig => zeig%next
enddo  ! zeig (cohorts)  

! Calculation of Water Supply
frtrel_1 = 1.

select case (flag_wurz) 
case (0)
  if (nroot_max .gt. 5) then 
     nroot3 = 5
  else
     nroot3 = nroot_max
  endif
case default
  nroot3=nroot_max
end select  




! layers with seedlings   
do j = 1,nroot3
   ! determination of resisctance coefficient
   
   select case (flag_wred)
   case(1)
      hred = fred1(j)
   case(2)
      hred = fred2(j)
   case(3)
      hred = fred3(j)
   case(4)
      hred = fred4(j)
   case(5)
      hred = 1.
   case(6)
      hred = 0.5
   case(7)
      hred = 0.25
   case(8)  
      hred = fred6(j)
   case(10)
      hred = fred7(j)
   case(11)
      hred = fred11(j)
   end select
   
   wat_res(j) = hred

   if (temps(j) .gt. -0.3)  then
      wat_at = max(wats(j) - wilt_p(j), 0.)     ! total available water per layer
      wat_ar = hred * wat_at                    ! total available water per layer with uptake resistance
	  hupt = 0.
   else
      wat_ar = 0.                              ! frost
      wat_at = 0.
      hupt   = 0.
   endif



! Distribution of Water Supply into the Cohorts
! Distribution of Fine Roots
  
   zeig => pt%first
   i = 1   ! cohort index
   do while (associated(zeig))
     if (zeig%coh%species .ne. nspec_tree+2) then     ! not for mistletoe
      frtrel = zeig%coh%frtrelc(j)
     	
     wat_ava = frtrel * wat_ar        ! available water per tree cohort and layer 

    	
	 if (wat_ava .ge. tr_dem(i)) then
         hupt_c    = tr_dem(i)
         tr_dem(i) = 0.
     else
         hupt_c    = wat_ava
         tr_dem(i) = tr_dem(i) - wat_ava
     endif
	     
     select case (flag_dis)    !xylem clogger disturbance 
     case (1,2)
     hupt_c = hupt_c * xylem_dis
     end select
	 
     xwatupt(i,j) = hupt_c   ! water uptake per cohorte and layer
     zeig%coh%supply = zeig%coh%supply + hupt_c
     if (zeig%coh%supply .lt.0.) then
     continue
     endif
     hupt = hupt + hupt_c
  

	 i = i + 1
    end if              ! exclusion of mistletoe
     zeig => zeig%next
   enddo  ! zeig (cohorts) 

   wupt_r(j) = hupt 

enddo    ! j

! layers without seedlings   

if (totfrt_p.gt.(seedlfrt+zero)) then
 totfrt_2 = 1./(totfrt_p-seedlfrt)

  do j = nroot3+1, nroot_max
   ! determination of resisctance coefficient
   select case (flag_wred)
   case(1)
      hred = fred1(j)
   case(2)
      hred = fred2(j)
   case(3)
      hred = fred3(j)
   case(4)
      hred = fred4(j)
   case(5)
      hred = 1.
   case(6)
      hred = 0.5
   case(7)
      hred = 0.25
   end select
   
   wat_res(j) = hred
   
   if (temps(j) .gt. -0.3)  then
      wat_at = max(wats(j) - wilt_p(j), 0.)     ! total available water per layer
      wat_ar = hred * wat_at                    ! total available water per layer with uptake resistance
      hupt   = 0.
   else
      wat_ar = 0.
   endif
   
   zeig => pt%first
   i = 1   ! cohort index
   do while (associated(zeig))
     frtrel = zeig%coh%frtrelc(j)		
     wat_ava = frtrel * wat_ar        ! available water per tree cohort and layer 
     if (wat_ava .ge. tr_dem(i)) then
         hupt_c    = tr_dem(i)
         tr_dem(i) = 0.
     else
         hupt_c    = wat_ava
         tr_dem(i) = tr_dem(i) - wat_ava
     endif
       
     select case (flag_dis)    !xylem clogger disturbance 
     case (1,2)
     hupt_c = hupt_c * xylem_dis
     end select
	 
	 xwatupt(i,j) = hupt_c
     zeig%coh%supply = zeig%coh%supply + hupt_c
     hupt = hupt + hupt_c
     i = i + 1
     zeig => zeig%next
   enddo  ! zeig (cohorts) 

   wupt_r(j) = hupt 
enddo     ! j
endif
END	subroutine upt_wat

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

SUBROUTINE upt_wat1

! Water uptake by roots
! 2. Version

use data_simul
use data_evapo
use data_soil
use data_stand 
use data_par

implicit none

real, dimension(1:anz_coh) :: tr_dem,frt_rel   ! help arrays for cohorts
real wat_ava, hdem, frtrel, frtrel_1, hupt, hupt_c, totfrt_2
real wat_at    ! total available water per layer
real wat_ar    ! total available water per layer with uptake resistance
real hred      ! resistance coefficient
real, external :: fred1, fred2, fred3, fred4, fred5, fred6, fred7, fred11

integer i, ianz, j, nroot3

ianz = anz_coh
tr_dem=0
trans_dem = (pet-aev_i) * alfm * (1. - exp(-gp_can/gpmax))  ! pet NOT reduced by ground evaporation
if (trans_dem .lt. 0.) trans_dem = 0.

! potential transpiration demand of each cohort
if (gp_can .gt. zero) then
   hdem   = trans_dem / gp_can
else
   hdem= 0.
endif

! Estimation of transpiration demand of tree cohorts and total fine root mass
! in layers with and without seedlings
zeig => pt%first
i = 1
do while (associated(zeig))
  select case (flag_eva)
  case (0, 1, 3)
    zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem  
  case (2)
    zeig%coh%demand = zeig%coh%demand - zeig%coh%aev_i
  end select
  tr_dem(i) = zeig%coh%demand
  i = i + 1
  zeig => zeig%next
enddo  ! zeig (cohorts)  

! uptake controlled by share of roots 
frtrel_1 = 1.

! layers with seedlings   
do j = 1,nroot_max
   ! determination of resisctance coefficient
   select case (flag_wred)
   case(1)
      hred = fred1(j)
   case(2)
      hred = fred2(j)
   case(3)
      hred = fred3(j)
   case(4)
      hred = fred4(j)
   case(5)
      hred = 1.
   case(6)
      hred = 0.5
   case(7)
      hred = 0.25
   case(8)  ! BKL, ArcEGMO
      hred = fred6(j)
   case(10)
      hred = fred7(j)
   end select

   wat_at = max(wats(j) - wilt_p(j), 0.)     ! total available water per layer
   wat_ar = hred * wat_at                    ! total available water per layer with uptake resistance
   hupt   = 0.
   
   zeig => pt%first
   i = 1   ! cohort index
   do while (associated(zeig))
     
     frtrel = zeig%coh%frtrel(j) * zeig%coh%x_frt * zeig%coh%ntreea * totfrt_1    
     wat_ava = frtrel * wat_ar        ! available water per tree cohort and layer 
     if (wat_ava .ge. tr_dem(i)) then
         hupt_c    = tr_dem(i)
         tr_dem(i) = 0.
     else
         hupt_c    = wat_ava
         tr_dem(i) = tr_dem(i) - wat_ava
     endif
	     
     select case (flag_dis)    !xylem clogger disturbance 
     case (1,2)
     hupt_c = hupt_c * xylem_dis
     end select
	 
     xwatupt(i,j) = hupt_c
     zeig%coh%supply = zeig%coh%supply + hupt_c
     hupt = hupt + hupt_c
     i = i + 1
     zeig => zeig%next
   enddo  ! zeig (cohorts) 

   wupt_r(j) = hupt 
enddo    ! j
END	subroutine upt_wat1

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

real FUNCTION fred1(j)

! Function for calculating uptake resistance
! from CHEN (1993)
! empirical relation between soil water content and resistance
! fred1=1 if (field_cap - 10%*field_cap) <= wats <= (field_cap + 10%*field_cap)

use data_par
use data_soil

implicit none
real hf, f09, f11, wc, diff
integer j

f09 = 0.9 * field_cap(j)
f11 = 1.1 * field_cap(j)
wc  = wats(j)

if (wc .lt. wilt_p(j)) then
   hf = 0.
else if (wc .lt. f09) then
   diff = f09-wilt_p(j)
   if (diff .lt. zero) diff = 0.001
   hf = 1. - (f09-wc) / diff
else if (wc .gt. f11) then
   diff = pv(j)-f11
   if (diff .lt. zero) diff = 0.001
   hf = 0.3 + 0.7 * (pv(j)-wc) / diff
   if (hf .lt. zero) hf = 0.001
else
   hf = 1.
endif
fred1 = hf
END	function fred1

!**************************************************************
 
real FUNCTION fred2(j)

! Function for calculating uptake resistance
! from Aber and Federer (f=0.04 fuer Wasser in cm)
! only 40% of total available water are plant available per day

implicit none

integer j

fred2 = 0.05 

END	function fred2

!**************************************************************
 
real FUNCTION fred3(j)

! Function for calculating uptake resistance
! from CHEN (1993)
! modified to a profile defined in fred 
! fred3 may be described by a function (old version): 
!     fred3 = 0.0004*j*j - 0.0107*j + 0.0735
! this case: set from a root profile, defined by input of root_fr

use data_par
use data_soil

implicit none
real hf, f09, f11, wc, diff
! hf is a reduction factor in dependence on water content
real fred(15)

integer j

! uptake reduction depending on water content
f09 = 0.9 * field_cap(j)
f11 = 1.1 * field_cap(j)
wc  = wats(j)

if (wc .lt. wilt_p(j)) then
   hf = 0.
else if (wc .lt. f09) then
   diff = f09-wilt_p(j)
   if (diff .lt. zero) diff = 0.001
   hf = 1. - (f09-wc) / diff
else if (wc .gt. f11) then
   diff = pv(j)-f11
   if (diff .lt. zero) diff = 0.001
   hf = 0.3 + 0.7 * (pv(j)-wc) / diff
   if (hf .lt. zero) hf = 0.001
else
   hf = 1.
endif

fred3 = root_fr(j) * hf 

END	function fred3

!**************************************************************
 
real FUNCTION fred4(j)

! Function for calculating uptake resistance
! modified to a profile defined in fred 
! profile at Beerenbusch

use data_soil

implicit none
real fred(15)

integer j

fred = (/ 0.0, 0.03, 0.03, 0.02, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01,    &
              0.01, 0.01, 0.01, 0.01, 0.01 /)  ! fred fuer Beerenbusch
 
fred4 = fred(j) 

END	function fred4

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

real FUNCTION fred6(j)

! Function for calculating uptake resistance
! from Kloecking (2006) simular to fred1
! empirical relation between soil water content and resistance
! fred6=1 if field_cap  <= wats <= (field_cap + 10%*field_cap)

use data_soil

implicit none
real hf, f09, f11, wc
integer j

f09 = field_cap(j)
f11 = 1.1 * field_cap(j)
wc  = wats(j)

if (wc .le. wilt_p(j)) then
   hf = 0.
else if (wc .lt. f09) then
   hf = 0.1 + (0.9 *(wc-wilt_p(j)) / (f09-wilt_p(j)))
else if (wc .gt. f11) then
   hf = 0.3 + 0.7 * (pv(j)-wc) / (pv(j)-f11)
   if (hf .lt. 0.) hf = 0.001
else
   hf = 1.
endif

fred6 = hf

END	function fred6

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

real FUNCTION fred7(j)

! Function for calculating uptake resistance
! from CHEN (1993)
! empirical relation between soil water content and resistance
! fred1=1 if (field_cap - 10%*field_cap) <= wats <= (field_cap + 10%*field_cap)

use data_par
use data_soil

implicit none
real hf, f09, f11, wc, diff
integer j

f09 = 0.9 * field_cap(j)
f11 = 1.1 * field_cap(j)
wc  = wats(j)

if (wc .lt. wilt_p(j)) then
   hf = 0.
else if (wc .lt. f09) then
   diff = f09-wilt_p(j)
   if (diff .lt. zero) diff = 0.001
   hf = exp(-5.*(f09-wc) / diff)
else if (wc .gt. f11) then
   diff = pv(j)-f11
   if (diff .lt. zero) diff = 0.001
   hf = 0.3 + 0.7 * (pv(j)-wc) / diff
   if (hf .lt. zero) hf = 0.001
else
   hf = 1.
endif

fred7 = hf

END	function fred7

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

real FUNCTION fred11(j)

! Function for calculating uptake resistance, especially adapted for Mistletoe disturbance
! function after van Wijk, 2000

use data_par
use data_soil
implicit none
real hf, S, f11, wc, diff
integer j
f11 = 1.1 * field_cap(j)
wc  = wats(j)
if (wc .lt. wilt_p(j)) then
   hf = 0.
else if (wc .lt. field_cap(j)) then
   S=(field_cap(j)-wc)/(field_cap(j)-wilt_p(j))
   hf = exp(-30*S)                                !30 = strong reduction in water avail.
else if (wc .gt. f11) then
   diff = pv(j)-f11
   if (diff .lt. zero) diff = 0.001
   hf = 0.3 + 0.7 * (pv(j)-wc) / diff
   if (hf .lt. zero) hf = 0.001
else
   hf = 1.
endif
fred11 = hf
END function fred11

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

 





 
 SUBROUTINE take_wat(eva_dem, psi)   
      
! Estimation of water taking out for uncovered soil 
use data_soil
use data_simul

implicit none

!input:
real :: eva_dem	   ! evaporation demand
real :: psi        ! covering

integer i, ii, j, ntag  ! max. layer of taking out
real, allocatable, dimension(:) :: gj
real, external :: b_r, funcov
real diff, gj_j, depth_j, depth_n, rij, rmax, rr, rs, sr

allocate (gj(nlay))

do i=1,nlay
   wupt_ev(i)=0.0
   gj(i)=0.0 
enddo  
      
ntag    = 0
rmax    = 0.0
depth_n = depth(n_ev_d)
	
do i=1,n_ev_d
  rij = 0.0
  rr  = depth_n/depth(i)
  sr  = 0.0
  rs  = 0.0
	    
  do j=1,i
!         depth for uncovered take out
     depth_j = depth(j)
     gj(j)   = FUNCOV(w_ev_d, rs*rr, rr*depth_j)
     rs      = depth_j
     sr      = sr + gj(j)
  enddo		! i
	  
  if (sr.gt.1.E-7) then
     sr = 1.0/sr
	    
     do j=1,i
!   	   water take out
!         (psi = 1.-psi) no soil evaporation in case of total covering 
!          and maximal evaporation for uncovered soil
        gj_j = -B_R(wats(j), field_cap(j), wilt_p(j))    &
               * eva_dem * (1.-psi) * gj(j) * sr
        gj_j = max(gj_j,0.0)
        gj(j)= gj_j
        rij  = rij + gj_j
     enddo		! i

     if (rij .gt. rmax) then
         rmax = rij
         ntag = i

         do ii=1,ntag
            wupt_ev(ii) = gj(ii)
         enddo	

     endif	! rij 
   endif	! sr	  
enddo	! n_ev_d

! balance
do i=1,nlay
   diff = wats(i) - wilt_p(i)
   if (wupt_ev(i) .gt. diff) then
       wupt_ev(i) = diff
   endif
enddo ! nlay

deallocate (gj)

END  subroutine  take_wat

!*******************************************************************************
   
real FUNCTION B_R(water, f_cap, wilting)

! Reduction function for water taking out (uncovered soil)

implicit none

!input:
real :: water	    ! water storage
real :: f_cap	    ! field capacity
real :: wilting	    ! wilting point

b_r  = 1.0

if (water .lt. f_cap) B_R = max((water-wilting)/(f_cap-wilting), 0.0)

END  function B_R

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

real FUNCTION funcov(wt_d, a, bb)

! take out density function for uncovered soil

implicit none

!input:
real :: wt_d    ! depth of water taking out by evaporation (cm)
real :: a, bb   ! relative upper and lower depth of actual layer
real fk, wt_5, b

      fk   = .455218234
      wt_5   = 0.05 * wt_d
      b      = min(bb,wt_d)
      funcov = (- b + a + 1.05*wt_d*log((b+wt_5)/(a+wt_5)))*fk/wt_d
      
END function funcov

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

real FUNCTION wat_new(wat_us, wat_in, ilayer)
!      FUNCTION WIEN(WIA,NIST,ALAM,DTI,TT,DICK)

! Estimation of additional water after infiltration and percolation

use data_par
use data_soil

implicit none

! input:
real :: wat_us	   ! water content in relation to field capacity
real :: wat_in	   ! water infiltration into actual layer
integer :: ilayer  ! number of actual layer
real dti	 !time step
real awi, b1, b2, la, hsqr, exphelp

dti = 1.
fakt = 0.4

if (fakt .ge. 0.0) then				! percolation?
   la = 100.0 * fakt * dti * wlam(ilayer)/thick(ilayer)**2
   if (wat_us .le. zero) then		! water near zero?
      if (wat_in .le. zero) then	! infiltrated water near zero?
         wat_new = wat_us + wat_in
      else
         if (wat_us+wat_in .gt. zero) then
            exphelp = sqrt(la*wat_in) * (1 + wat_us/wat_in)*1
            if (exphelp .le.10.) then     ! avoid underflow   
               b1 = -exp(-2. * exphelp)
            else
               b1 = 0.
            endif
            wat_new = sqrt(wat_in/la) * (1+b1)/(1-b1)
		 else
		    wat_new = wat_us + wat_in
		 endif
	  endif	! wat_in

   else 
        if (wat_in .lt. 0.) then
            awi = abs(wat_in)
            b1 = atan(wat_us/sqrt(awi/la)) / sqrt(la * awi)
			if (b1 .gt. 1) then
			   b2 = sqrt (awi * la)
			   b1 = sin(b2) / cos(b2)
			   b2 = sqrt(awi / la)
			   wat_new = b2 * (wat_us - b2*b1) / (b2 + wat_us*b1)
			else
			   wat_new = wat_in * (1-b1)
			endif	! b1
        else 
             if (wat_in .gt. 0.) then
                b1 = sqrt(wat_in / la)
                hsqr = sqrt(la*wat_in)
                if (hsqr .lt. 10.) then
                   b2 = (wat_us - b1) * exp(-2.* hsqr) / (wat_us + b1)
                   if (b2 .ge. 1.0) then
                       b2 = 0.99999
                   endif
                else
                   b2 = 0.
                endif
                wat_new = b1 * (1.+b2) / (1.-b2)
             else
                wat_new = wat_us / (1. + la*wat_us)
             endif	
       endif	! wat_in		   
   endif	! wat_us
else
    wat_new = wat_us
endif		! fakt

END  function wat_new
 
!******************************************************************************

SUBROUTINE bucket(bucksize1, bucksize2, buckdepth)

! calculation of bucket size (1m; without humus layer)

use data_soil

implicit none

real bucksize1,   &    ! bucket size of 1 m depth  (nFK)
     bucksize2,   &    ! bucket size of rooting zone
     buckdepth, diff
integer j

bucksize1 = 0.
bucksize2 = 0.
buckdepth = 0.
do j=2,nlay
   if ((depth(j)-depth(1)) .lt. 100.) then
      bucksize1  = bucksize1 + wats(j) - wilt_p(j)
	  buckdepth = depth(j) - depth(1)
   else
      diff      = 100. - buckdepth
	  bucksize1  = bucksize1 + (wats(j) - wilt_p(j))*diff/thick(j)
      buckdepth = 100.
	  exit
   endif
enddo

do j=2,nroot_max
      bucksize2  = bucksize2 + wats(j) - wilt_p(j)
enddo

END  subroutine bucket

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

SUBROUTINE snowpack(snow_sm, p_inf, pev)

! properties of snow 
! calculation of soil surface temperature under snow pack 

use data_climate
use data_evapo
use data_inter
use data_par
use data_simul
use data_soil
use data_soil_t

implicit none

real p_inf		   ! infiltrated water
real snow_sm
real pev
real airtemp_sm    ! melting temperature 
real snow_old      ! old snow pack
real tc_snow       ! thermal conductivity of snow  J/cm/s/K
real thick_snow    ! thickness of snow
real dens_snow     ! density of snow
real:: dens_sn_new  = 0.1 ! density of fresh snow
real fakta

snow_old = snow

!substract evaporation of snowcover from snow in both cases
if (airtemp .lt. temp_snow) then        ! frost conditions
    snow     = snow + prec_stand	    ! precipitation as snow
    snow_sm  = 0.0          	        ! no snow melting
    p_inf    = 0.0                      ! no infiltrated precipitation
    pev	     = max((pev_s - aev_i), 0.) ! interc. evapor. reduces soil evapor.

else

    airtemp_sm = max(airtemp, 0.)
    snow_sm    = airtemp*(0.45+0.2*airtemp)  ! snow melting
    snow_sm    = MIN(snow_sm, snow)
    snow       = snow - snow_sm
    p_inf      = prec_stand + snow_sm	  ! infiltrated precipitation	
    pev	       = max((pev_s - aev_i), 0.) ! interc. evapor. reduces soil evapor.

end if		! airtemp

if (snow .ge. zero) then
   snow_day = snow_day + 1
   days_snow = days_snow + 1
   if (pev .le. zero) then
      pev = 0.	
   else
    ! snow sublimation 
      aev_s = max(min(snow, pev), 0.)   
      snow  = snow - aev_s
      pev   = pev - aev_s
   endif

   ! soil surface temperature under snow pack
   ! snow hight = 0.2598 * water equivalent + 8.6851; adjustment from measurement values (see Bodentemperatur.xls)

   if (snow .ge. 0.05) then
      
      dens_snow = dens_sn_new + snow_day*0.025
      dens_snow = MIN(dens_snow, 1.)
      dens_snow = 0.5*(dens_sn_new*prec_stand + dens_snow*snow_old)/snow
      dens_snow = MIN(dens_snow, 1.)
      tc_snow   = 0.7938*EXP(3.808*dens_snow)*0.001  ! thermal conductivity of snow  J/cm/s/K
      thick_snow = snow / dens_snow 
      fakta      = tc_snow * 86400. * (thick(1)/2.) / (t_cond(1) * thick_snow)   ! s --> day
      temps_surf  = (0.5*temps(1) + fakta*airtemp) / (1. + fakta)   ! CoupModel (Jansson, 2001)
   endif   					
else
   snow_day = 0
endif
END  subroutine snowpack

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

SUBROUTINE soil_stress 

! Calculation of the stress factors

use data_soil
use data_species
use data_stand
use data_par

implicit none
integer :: i, k

real :: m_1, m_2, n_1, n_2
real :: wratio, wafpo
real, dimension (1:4) ::  allstress, xvar, yvar

!temperature stress
do i=1,nlay
 do k=1,nspecies
    if (temps(i) .ge. spar(k)%tbase) then
        svar(k)%tstress(i) = sin((pi/2)*(temps(i)-spar(k)%tbase)/(spar(k)%topt-spar(k)%tbase))
   	else
	    svar(k)%tstress(i) = 0.
    endif

 !soil strength
    wratio=0.    
    if (dens(i) .le. BDopt(i)) then
		svar(k)%BDstr(i) = 1
		svar(k)%BDstr(i) = 1
	elseif (dens(i) .ge. svar(k)%BDmax(i)) then
		svar(k)%BDstr(i) = 0
	else
		svar(k)%BDstr(i) = (svar(k)%BDmax(i)-dens(i))/(svar(k)%BDmax(i)-BDopt(i))
	endif

	if (watvol(i) .lt. wilt_p_v(i)) then
	wratio = 0.
	elseif (watvol(i) .gt. f_cap_v(i)) then
	wratio = 1.
	else
	wratio = (watvol(i)-wilt_p_v(i))/(f_cap_v(i)-wilt_p_v(i))
	endif
	
	svar(k)%sstr(i)=svar(k)%BDstr(i)*sin(1.57*wratio)
 
 !aeration
    wafpo=watvol(i)/pv_v(i)
	if (wafpo .ge. svar(k)%porcrit(i)) then
		svar(k)%airstr(i) = (1.-wafpo)/(1.-svar(k)%porcrit(i))
	else
		svar(k)%airstr(i) = 1.
	endif
	
    if (svar(k)%airstr(i) .lt. 0.) svar(k)%airstr(i) = 0.
	
 !soil acidity
    xvar=(/spar(k)%ph_min, spar(k)%ph_opt_min, spar(k)%ph_opt_max, spar(k)%ph_max/)
    yvar=(/0,1,1,0/)
    m_1=(yvar(1)-yvar(2))/(xvar(1)-xvar(2))
    n_1=yvar(2)-m_1*xvar(2)
    m_2=(yvar(3)-yvar(4))/(xvar(3)-xvar(4))
    n_2=yvar(4)-m_2*xvar(4)

	if (phv(i) .gt.  spar(k)%ph_opt_max .and. phv(i) .le. spar(k)%ph_max ) then
		svar(k)%phstr(i)=m_2*phv(i)+n_2
		elseif (phv(i) .lt. spar(k)%ph_opt_min .and. phv(i) .ge. spar(k)%ph_min ) then
			svar(k)%phstr(i)=m_1*phv(i)+n_1
			elseif (phv(i) .gt. spar(k)%ph_max .or. phv(i) .lt. spar(k)%ph_min) then
			svar(k)%phstr(i)=0.
				else
				svar(k)%phstr(i)=1.
	endif

 ! total stress (Rstress) is taken as the largest of the four
    
	allstress(1)=svar(k)%tstress(i)
	allstress(2)=svar(k)%sstr(i)
	allstress(3)=svar(k)%airstr(i)
	allstress(4)=svar(k)%phstr(i)
	
	svar(k)%Rstress(i)= minval(allstress)
	svar(k)%Smean(i)=svar(k)%Rstress(i)+svar(k)%Smean(i)
	
 enddo
enddo

END subroutine soil_stress

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

SUBROUTINE hum_add(xfcap, xwiltp, xpv) 
! Soil parameter according to [Kuntze et al., Bodenkunde, 1994], S. 172 

use data_simul
use data_soil
use data_soil_cn

implicit none
integer :: i, k
real    :: fcapi, clayvi, siltvi, humvi, humvi2, wiltpi, pvi, nfki, hcbc 
real, dimension(nlay):: xfcap, xwiltp, xpv   ! output of addition mm/dm

    xfcap(1)  = 0.0
    xwiltp(1) = 0.0
    xpv(1)    = 0.0

do i = 1, nlay
    fcapi  = 0.
    wiltpi = 0.
    pvi    = 0.
    clayvi = clayv(i)
    humvi  = humusv(i)*100.
    humvi2 = humusv(i)*humusv(i)
    if (humvi .lt. 15.) then
        if (clayvi .le. 0.05) then
            wiltpi = 0.0609 * humvi2 + 0.33 * humvi 
            pvi    = 0.0436 * humvi2 + 0.631 * humvi
            nfki   = -0.0009 * humvi2 + 1.171 * humvi
            fcapi  = nfki + wiltpi
        
        else if (clayvi .le. 0.12) then
            wiltpi = 0.0357 * humvi2 + 0.0762 * humvi 
            pvi    = 0.0441 * humvi2 + 0.5455 * humvi
            nfki   = 0.0252 * humvi2 + 0.7462 * humvi
            fcapi  = nfki + wiltpi
            
        else if (clayvi .le. 0.17) then
            wiltpi = 0.0374 * humvi2 - 0.1777 * humvi 
            pvi    = 0.0552 * humvi2 + 0.2936 * humvi
            nfki   = 0.0324 * humvi2 + 0.6243 * humvi
            fcapi  = nfki + wiltpi
        
        else if (clayvi .le. 0.35) then
            wiltpi = 0.0179 * humvi2 - 0.0385 * humvi 
            pvi    = 0.0681 * humvi2 + 0.0768 * humvi
            nfki   = 0.0373 * humvi2 + 0.3617 * humvi
            fcapi  = nfki + wiltpi
       
        else if (clayvi .le. 0.65) then
            wiltpi = 0.0039 * humvi2 + 0.0254 * humvi 
            pvi    = 0.0613 * humvi2 + 0.0947 * humvi
            nfki   = 0.0338 * humvi2 + 0.0904 * humvi
            fcapi  = nfki + wiltpi
        
        else 
            wiltpi = 0.0 
            pvi    = 0.0613 * humvi2 + 0.0947 * humvi
            nfki   = 0.0104 * humvi2 + 0.2853 * humvi
            fcapi  = nfki + wiltpi
        
        endif     
    else   ! humvi > 15
       ! organic soils
        continue
    endif   ! humvi

    xfcap(i)  = fcapi
    xwiltp(i) = wiltpi
    xpv(i)    = pvi
enddo
    
if (flag_bc .gt. 0) then
    do i = 1, nlay
        if (C_bc(i) .gt. 0.) then
            fcapi  = f_cap_v(i)
            clayvi = clayv(i)
            siltvi = siltv(i)
            humvi  = humusv(i)*100.
            hcbc   = C_bc(i)*100.*100. / (cpart_bc(y_bc_n) * dmass(i))
            if ((clayvi .le. 0.17) .and. (siltvi .le. 0.5)) then     ! sand
                fcapi  = 0.0619 * hcbc
                wiltpi = 0.0375 * hcbc
                nfki = 7.0   
            elseif ((clayvi .le. 0.45) .and. (siltvi .gt. 0.17)) then   ! loam
                fcapi  = 0.015 * hcbc
                wiltpi = 0.0157 * hcbc
                nfki = 10.
            else      ! clay
                fcapi  = -0.0109 * hcbc
                wiltpi = -0.0318 * hcbc
                nfki = 16.
            endif
            xfcap(i)  = xfcap(i) + fcapi
            xwiltp(i) = xwiltp(i) + wiltpi
       endif       
    
    enddo
endif

END subroutine hum_add

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

SUBROUTINE bc_appl 

! application of biochar

use data_out
use data_simul
use data_soil
use data_soil_cn

implicit none

character :: text
integer :: ios, inunit, j
logical :: ex
real    :: hcbc 

    call testfile(valfile(ip),ex)
    IF (ex .eqv. .true.) then
      inunit = getunit()
      ios=0
      open(inunit,file=valfile(ip),iostat=ios,status='old',action='read')
      if (.not.flag_mult8910) then
          print *,'***** Reading application values of biochar from file ',valfile(ip),'...'
          write (unit_err, *) 'Application values of biochar from file   ',trim(valfile(ip))
      endif

      do
         read(inunit,*) text
         IF(text .ne. '!')then
            backspace(inunit)
            exit
         endif
      enddo
      
      read (inunit,*,iostat=ios) n_appl_bc
      allocate (C_bc_appl(n_appl_bc))
      allocate (N_bc_appl(n_appl_bc))
      allocate (bc_appl_lay(n_appl_bc))
      allocate (cnv_bc(n_appl_bc))
      allocate (dens_bc(n_appl_bc))
      allocate (cpart_bc(n_appl_bc))
      allocate (y_bc(0 : n_appl_bc + 1))
      y_bc = 0
      C_bc_appl = 0.
      N_bc_appl = 0.
      do j = 1, n_appl_bc
          read (inunit,*,iostat=ios) y_bc(j), cpart_bc(j), cnv_bc(j), dens_bc(j)
          read (inunit,*,iostat=ios) bc_appl_lay(j), C_bc_appl(j)      
      enddo
    endif   ! ex

END subroutine bc_appl

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