Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • foresee/4C
  • gutsch/4C
2 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 8559 additions and 0 deletions
source_code/version_2.3_windows/simul.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* - Simulation control: SIM_CONTROL *!
!* SIMULATION_4C *!
!* *!
!* 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 sim_control
use data_climate
use data_simul
use data_site
use data_out
implicit none
integer run_nr, ipp, irl, icl, i
character a
character(8) actdate
character(10) acttime, helpsim, text1, text2
real time1, time2, time3
unit_err=getunit()
if(flag_multi.eq.5) dirout = './'
open(unit_err,file=trim(dirout)//trim(site_name(1))//'_error.log',status='replace', position='append')
unit_trace=getunit()
open(unit_trace,file=trim(dirout)//trim(site_name(1))//'_trace.log',status='replace', position='append')
write (unit_trace, *) ' Trace of calls - subroutines of 4C '
write (unit_trace, *)
write (unit_trace, *) 'iday time_cur subroutine '
write (unit_trace, '(I4,I10,A)') iday, time_cur, ' sim_control'
! check daily output
if (year > 5 .and. flag_dayout .ge. 1) then
write(*,*) ' Warning: Your choice of daily output is ON with a simulation time of'
write(*,'(I6,A,I8,A)') year,' years. This option will create ',365*year,' data records per file!'
write(*,'(A)',advance='no') ' Do you really want do use daily output (y/n)? '
read *,a
IF (a .eq. 'n' .or. a .eq. 'N') then
flag_dayout = 0
ENDIF
ENDIF
! open file ycomp (yearly compressed output (multi run))
IF (time_out .ne. -2) call prep_out_comp
!call epsilon
IF(flag_multi.eq.1) THEN
run_nr = site_nr
ELSE IF (flag_multi.eq.5) THEN
run_nr = 1
ELSE
run_nr = repeat_number
ENDIF
call date_and_time(actdate, acttime)
write (unit_err, *)
time3 = 0.
if (.not.flag_mult910) then
nrreal = 1
nrclim = 1
endif
do icl = 1, nrclim ! climate scenarios
iclim = icl
DO ipp = 1, run_nr ! sites
ip = ipp
if (flag_trace) write (unit_trace, '(I4,I10,A,I3)') iday, time_cur, ' sim_control ip=',ip
do irl = 1, nrreal ! realization f climate scenarios
if (flag_mult910) then
climfile(ip) = climszenfile(ip, icl, irl)
site_name (ip) = trim(site_name1)//'_'//trim(sitenum(ip))
write (helpsim,'(I10)') icl
read (helpsim,*) text1
write (helpsim,'(I10)') irl
read (helpsim,*) text2
site_name (ip) = trim(site_name (ip))//'_'//trim(text1)//'_'//trim(text2)
write (unit_err, *)
write (unit_err, '(A,3I5,2X, A50)')'* ip, cli-scenario, real., site: ', ip, icl, irl, site_name(ip)
write (unit_err, '(A,A)') 'Climate file: ', trim(climfile(ip))
else
write (unit_err, *)
write (unit_err, '(A10,I5,2X, A50)') ' ip/site ', ip, trim(site_name(ip))
site_name1 = trim(site_name(ip))
endif
call CPU_time (time1)
if(ip.ne.0) then
CALL sim_ini
CALL prepare_site
if (flag_multi.eq.5) then
! call m4c_simenv_in
unit_comp2 = 6 ! standard output
end if
if(flag_end.gt.0) then
select case (flag_end)
case (1)
print*,ip, ' stop in prepare_stand (see error.log)'
case (2)
print*, ip, 'stop in prepare_stand, stand ', &
'identificator not found in prepare_stand'
case (3)
print*,ip, 'stop in canopy'
case (4)
print*,ip, 'stop in readsim, climate ID not found'
case (5)
print*,ip, ' stop in readsoil, soil ID not found ', adjustl(soilid(ip))
case (6)
write(*,'(A,I5)') ' >>>foresee message: stop in read_cli, no climate data for year ',time_b
call finish_simul
stop
case (7)
print*,ip, ' stop in readsoil, error during reading soil data ', adjustl(soilid(ip))
call finish_simul
stop
case default
print*,ip, 'flag_end = ', flag_end
end select
call finish_simul
flag_end = 0
else
IF(flag_multi==2) CALL fixclimscen
if (.not.flag_mult910) then
write (*,*)
write (*,*) '>>> Start FORESEE-Simulation site ', ipp
write (*,*)
endif
CALL simulation_4c
CALL finish_simul
endif
if (flag_mult910) then
call out_var_stat(1, irl)
else
if ((flag_multi .ne. 8) .and. (nvar .gt. 1)) call out_var_stat(3, 1)
endif
endif ! flag_end
call CPU_time (time2)
if (.not.flag_mult910) then
print *, ' run time for simulation ',ip, time2-time1, ' sec'
endif
write (unit_err, *) ' run time for simulation ',ip, time2-time1, ' sec'
time3 = time3 + (time2-time1)
enddo ! irl
if (flag_mult910) call out_var_stat(2, -99)
write (unit_err, *)
write (unit_err, *)
write (unit_err, *) '* * * * * New ip/site * * * * *'
ENDDO ! ip until site_nr (page number)
write (unit_err, *)
write (unit_err, *) '************ New climate scenario **********'
enddo ! icl
if (nvar .gt. 1) then
select case (flag_multi)
case (5, 9, 10)
continue
case (1)
continue
case default
call out_var_file
end select
endif
! comparison with measurements
if (flag_stat .gt. 0) CALL mess
call CPU_time (time1)
time3 = time3 + (time1-time2)
write (unit_err, *)
write (unit_err, *) ' total run time ', time3, ' sec'
CALL finish_all
PRINT *,'All simulations finished!'
END SUBROUTINE sim_control
!**************************************************************
SUBROUTINE simulation_4C
!*** Declaration part ***!
USE data_simul
USE data_species ! species specific parameters
USE data_site ! site specific data
USE data_climate ! climate data
USE data_soil
USE data_soil_cn
USE data_stand ! state variables of stand, cohort and cohort element
USE data_out
USE data_manag
USE data_plant
USE data_par
IMPLICIT NONE
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' simulation_4C'
! allocation of environmental variable fields
if(flag_wpm.ne.4 .and. flag_wpm.ne.5.and.flag_wpm.ne.6) then
! time loop
DO time = 1, year
iday = 1
! Update population variable for new year if population is changed through interventions
if (flag_standup .gt. 0 .or. flag_dis==1 .or. flag_dis==1) then
call stand_balance
call standup
flag_standup = 0
endif
CALL year_ini
! Calculate RedN from soil C/N
! read or create Redn for areal application
IF (time.EQ.1 .and. flag_redn) CALL RedN_ini
IF (flag_dis .eq. 1 .or. flag_dis .eq. 2) then
CALL dis_manag
endif
! simulation of processes with subannual resolution (fluxes and soil)
CALL stand_daily
if(flag_end.ge.1) exit ! exit do loop time
! compressed output of start values
IF (lcomp1) THEN
CALL out_comp(unit_comp1)
lcomp1 = .FALSE.
ENDIF
! cohort litter production
CALL senescence
! calculation of stand variables over all patches
CALL stand_balance
! calculation of soil variables for yearly output
CALL s_year
! calculation of fire variables for yearly output
CALL fire_year
! calculation of indices for yearly output
CALL t_indices(temp_mon)
! summation output
IF(flag_sum.eq.4) THEN
write(unit_sum,'(I5,9F11.3)') time_cur,photsum,npppotsum,nppsum,resosum,lightsum,nee,abslightsum,precsum, tempmean
photsum=0.;npppotsum=0.;nppsum=0.;resosum=0.;lightsum=0.;nee=0.;abslightsum=0.;precsum=0.
ENDIF
totsteminc = 0.
totsteminc_m3 = 0.
! cohort loop for change in crown dimensions, allocation and tree dimension calculations
zeig=>pt%first
DO
IF (.not.ASSOCIATED(zeig)) exit
IF (zeig%coh%height.ge.thr_height .and. zeig%coh%species.le.nspec_tree) then
! determine crown movement
CALL CROWN( zeig )
! allocate NPP to the various tree compartments
CALL PARTITION( zeig )
if(flag_end.ge.1) exit ! do loop
ENDIF
IF (zeig%coh%species.EQ.nspec_tree+1) then ! Ground vegetation
! allocate NPP to the various ground vegetation compartments
CALL PARTITION_SV( zeig )
ENDIF
IF (zeig%coh%species.eq.nspec_tree+2) then ! Mistletoe
CALL PARTITION_MI( zeig )
if(flag_end.ge.1) exit ! do loop
ENDIF
zeig=>zeig%next
END DO
if(flag_end.ge.1) exit ! exit do loop time
! calculation of annual mortality
IF(flag_mort.ge.1) CALL stand_mort
! annual growth of trees below thr_height, which are initialized by planting (not seeded!)
! at the beginning of the simulation or during management (shelter-wood)
if(flag_reg.ne.2.and.flag_sprout.eq.0) CALL growth_seed
CALL mort_seed
if(flag_sprout.eq.1) flag_sprout=0
IF(flag_mg==1) then
if(thin_year(act_thin_year)==time_cur) then
CALL management
act_thin_year = act_thin_year+1
end if
ELSE IF((flag_mg.ge.2 .or. flag_mg.eq.3 .or. flag_mg.eq.33.or. flag_mg.eq.9 .or. flag_mg.eq.10 .or. flag_mg.eq. 333).and.anz_spec.ne.0) THEN
CALL management
if(flag_wpm.ne.0) CALL timsort
ENDIF
! no assortment if wpm is not called
if(flag_mg.eq.0.and.anz_spec.ne.0) then
if(flag_wpm.ne.0) call timsort
end if
CALL litter
! input of dead biomass into soil compartments
CALL cn_inp
! if(flag_multi.eq.5) call m4c_simenv_out
! annual establishment for all species
IF (flag_reg.eq.1.or.flag_reg.eq.2.or.flag_reg.eq.3.or.flag_reg.eq.30) CALL stand_regen
! cumsteminc = cumsteminc + totsteminc
! planting of seedlings/saplings at the beginning of simulation
if(flag_reg.ge.9..and. flag_reg.lt.100. .and. time.eq.1) call planting
if(flag_reg.ge.9..and. flag_reg.lt.100. .and. flag_mg .eq.44) call planting
! Update stand variables if stand changed
if (flag_standup.gt.0 .or. anz_spec.eq.0) then
call stand_balance
! if (flag_standup .gt. 1) call root_distr ! wird generell in year_ini berechnet
endif
cumsteminc = cumsteminc + totsteminc
! yearly output
IF (time_out .gt. 0) THEN
IF (mod(time,time_out) .eq. 0) then
CALL outyear (1)
CALL outyear (2)
endif
ENDIF
! store of output variables (multi run 4, 8, 9)
IF (nvar .gt. 1) CALL outstore
! RedN calculation
if ((flag_limi .eq. 10) .or. (flag_limi .eq. 15)) call RedN_calc
! CALL list_cohort
CALL del_cohort
if (.not.flag_mult910) PRINT *, ' * Year ', time, time_cur,' finished... '
END DO ! time
! calculation of stand variables over all patches at the end!
CALL stand_balance
!***** wpm ******
! check if management
if(flag_mg == 0) then
flag_wpm = 0
endif
if (flag_wpm == 1 .or. flag_wpm == 21 .or. flag_wpm == 11) call wpm
if (flag_wpm == 2) call sea
if (flag_wpm == 3) then
call wpm
call sea
end if
!*** * * * * * * * ****
else
call wpm
end if
if (flag_wpm .gt. 0) call out_wpm(1)
CALL out_comp(unit_comp2)
if(flag_end.eq.1) print*,ip, 'stop in partitio'
if(flag_end.eq.3) print*,ip, 'stop in calc_la in canopy: toplayer = 125 m'
flag_end = 0
if (.not.flag_mult910) PRINT *, ' * Simulation ',ip,' finished.'
END SUBROUTINE simulation_4C
!**************************************************************
source_code/version_2.3_windows/soil.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 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
!*******************************************************************************
source_code/version_2.3_windows/soil_cn.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 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 5C
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
!**************************************************************
source_code/version_2.3_windows/soil_cn_link.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* linking SOIL_C/N - Programs with forest module *!
!* *!
!* Author: F. Suckow *!
!* *!
!* contains: *!
!* S_CN_INI *!
!* S_CN_GENER *!
!* CN_INP *!
!* N_UPT(jlay) *!
!* READ_LITTER_INPUT *!
!* *!
!* 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 s_cn_ini
! Initialisation of soil data and parameters for C/N-module
use data_simul
use data_soil
use data_soil_cn
use data_species
use data_stand
implicit none
integer i, j
type (species_litter) :: sliti
real, external :: f_cnv, rmin_p, rnit_p
real :: xx, xcnv
! turnover biochar
k_bc = 0.00005
k_syn_bc =0.03
do j = 1, nlay
if (C_hum(j) .lt.0.) then
if (.not.flag_mult8910) call error_mess(time, 'missing value of C_hum set to 0.0 in layer ', real(j))
C_hum(j) = 0.0
endif
if (N_hum(j) .lt.0.) then
if (.not.flag_mult8910) call error_mess(time, 'missing value of N_hum set to 0.0 in layer ', real(j))
N_hum(j) = 0.0
endif
enddo
!!! zum Test ohne Primaersubstanz !!!
C_opm = 0.
N_opm = 0.
!!! zum Test ohne Primaersubstanz !!!
call s_cn_gener
! Convert concentration (mg/l) into contents (g/m2) per layer
NH4 = NH4 * 0.001 *wats
NO3 = NO3 * 0.001 *wats
if (flag_lit .eq. 0) then
! to even out balance for generated values ==> new values for C_ / N_hum
do j = 1, nlay
xx = C_hum(j)
if (N_hum(j) .gt. 1E-6) then
xcnv = f_cnv(C_hum(j), N_hum(j))
if (xx .gt. C_opm(j)) then
C_hum(j) = xx - C_opm(j)
N_hum(j) = C_hum(j) / xcnv
endif
endif
enddo
endif
! reduction of mineralization and nitrification depending on pH
do j=1,nlay
if (phv(j) .lt. 0) then
rmin_phv(j) = 1
rnit_phv(j) = 1
else
rmin_phv(j) = rmin_p(phv(j))
rnit_phv(j) = rnit_p(phv(j))
endif
cnv_opm(j) = f_cnv(C_opm(j), N_opm(j))
cnv_hum(j) = f_cnv(C_hum(j), N_hum(j))
enddo
call s_year ! calculate a year's values for start year as well
wats(1) = field_cap(1) ! ensuring consistency, in case novel calculation was done in s_year
! yearly cumulative quantities
N_min = 0.
N_lit = 0.
C_lit = 0.
C_accu = C_tot
Nleach_c = 0.
Nupt_c = 0.
Nupt_d = 0.
resps_c = 0.
END subroutine s_cn_ini
!**************************************************************
SUBROUTINE s_cn_gener
! Initialisation of soil data and parameters for C/N-module
use data_par
use data_simul
use data_soil
use data_soil_cn
use data_species
use data_stand
implicit none
integer i, j
real dbm_c ! C content of dead biomass
real dbm_frt ! C content of dead fine root biomass
real dbm_part ! part of dead biomass of previous years
real e_part ! part of dead biomass of one year
real t_day
real hconvd ! conversion factor kg/patchsize ==> g/m2
real hconvda ! conversion in C content and from tree to cohort
real, external :: f_cnv
type (species_litter) :: sliti
C_opm = 0.
N_opm = 0.
do i = 1, nspecies
if (i .eq. nspec_tree+2) then
continue
endif
sliti = slit(i)
sliti%species_name = spar(i)%species_name
if (flag_lit .eq. 0) then
sliti%C_opm_fol = 0.
sliti%C_opm_tb = 0.
sliti%C_opm_stem = 0.
sliti%C_opm_frt = 0.
sliti%C_opm_crt = 0.
endif
sliti%N_opm_fol = 0.
sliti%N_opm_tb = 0.
sliti%N_opm_stem = 0.
sliti%N_opm_frt = 0.
sliti%N_opm_crt = 0.
slit(i) = sliti
enddo
hconvd = 1000. / kpatchsize
if (flag_lit .eq. 0) then
zeig => pt%first
do
if (.not. associated(zeig)) exit
i = zeig%coh%species
if (i .ne. nspec_tree+2) then ! no litter initialisation for Mistletoe
sliti = slit(i)
hconvda = cpart * zeig%coh%ntreea
! consider decomposition rate, i.e. biomass of previous years
j = 1
t_day = 365.
dbm_part = 0.
do
! consider dependency on temp. and water
e_part = exp (-spar(i)%k_opm_fol * 0.2 * j * t_day)
dbm_part = dbm_part + e_part
if (e_part .gt. 0.001) then
j = j+1
else
exit
endif
enddo
select case (flag_dis)
case (1,2)
zeig%coh%litC_fol = spar(i)%psf * zeig%coh%x_fol * hconvda !(spar(i)%psf * zeig%coh%x_fol+zeig%coh%x_fol_loss) * hconvda, diese idee wieder weg ! Umrechnung in g/m2 erst in subr. litter
case (0)
zeig%coh%litC_fol = spar(i)%psf * zeig%coh%x_fol * hconvda ! conversion in g/m2 first into subr. litter
end select
zeig%coh%litN_fol = zeig%coh%litC_fol * (1.-spar(i)%reallo_fol) / spar(i)%cnr_fol
dbm_c = dbm_part * zeig%coh%litC_fol * hconvd
sliti%C_opm_fol = sliti%C_opm_fol + dbm_c
!dead fine root biomass of humus layer
! consider decomposition rate, i.e. biomass of previous years
j = 1
t_day = 365.
dbm_part = 0.
do
! consider dependency on temp. and water
e_part = exp (-spar(i)%k_opm_frt * 0.2 * j * t_day)
dbm_part = dbm_part + e_part
if (e_part .gt. 0.001) then
j = j+1
else
exit
endif
enddo
! change see foliage
select case (flag_dis)
case (1,2)
zeig%coh%litC_frt = spar(i)%psr * zeig%coh%x_frt * hconvda ! conversion in g/m2 first into subr. litter
case (0)
zeig%coh%litC_frt = spar(i)%psr * zeig%coh%x_frt * hconvda ! conversion in g/m2 first into subr. litter
end select
zeig%coh%litN_frt = zeig%coh%litC_frt * (1.-spar(i)%reallo_frt) / spar(i)%cnr_frt
dbm_c = dbm_part * zeig%coh%litC_frt * hconvd * (1.-spar(i)%reallo_frt) / spar(i)%cnr_frt
dbm_frt = dbm_c * zeig%coh%frtrel(1)
sliti%C_opm_frt(1) = sliti%C_opm_frt(1) + dbm_frt
sliti%N_opm_frt(1) = sliti%N_opm_frt(1) + dbm_frt
! Total fine root biomass must be distributed over all soil layers
do j = 2, nlay
dbm_frt = dbm_c * zeig%coh%frtrel(j)
sliti%C_opm_frt(j) = sliti%C_opm_frt(j) + dbm_frt
sliti%N_opm_frt(j) = sliti%N_opm_frt(j) + dbm_frt
enddo
slit(i) = sliti
endif ! (i .ne. nspec_tree+2)
zeig => zeig%next
enddo
endif
do i = 1, (nspec_tree+1) !exclusion of mistletoe
sliti = slit(i)
if (flag_lit .gt. 0) then
dbm_frt = sliti%C_opm_frt(1)
dbm_c = sliti%C_opm_crt(1)
do j = 1, nlay
sliti%C_opm_frt(j) = dbm_frt * root_fr(j)
sliti%N_opm_frt(j) = sliti%C_opm_frt(j) * (1.-spar(i)%reallo_frt) / spar(i)%cnr_frt
sliti%C_opm_crt(j) = dbm_c * root_fr(j)
enddo
endif
sliti%N_opm_fol = sliti%C_opm_fol * (1.-spar(i)%reallo_fol) / spar(i)%cnr_fol
! pools of dead biomass without stems
C_opm(1) = C_opm(1) + sliti%C_opm_fol + sliti%C_opm_tb + sliti%C_opm_frt(1) + sliti%C_opm_crt(1)
N_opm(1) = N_opm(1) + sliti%N_opm_fol + sliti%N_opm_tb + sliti%N_opm_frt(1) + sliti%N_opm_crt(1)
slit(i) = sliti
enddo
do j=2,nlay
C_opm(j) = SUM(slit%C_opm_frt(j)) ! + SUM(slit%C_opm_crt(j))
N_opm(j) = SUM(slit%N_opm_frt(j)) ! + SUM(slit%N_opm_crt(j))
if (C_opm(j) < 0.) then
continue
endif
enddo
! Total OPM of all species
END subroutine s_cn_gener
!**************************************************************
SUBROUTINE cn_inp
! Input of dead biomass (all fractions) into soil C- and N-pools
! call from simulation_4C
use data_simul
use data_par
use data_soil
use data_soil_cn
use data_species
use data_stand
implicit none
integer i, j
real hconvd, hf, hc, hfc, hfn, hfrtc, hfrtn, hfc1, Copm, Nopm, Clitf, Nlitf
type (species_litter) :: sliti
real, external :: f_cnv
Clitf = 0.
Nlitf = 0.
N_lit = 0.
C_lit = 0.
C_lit_fol = 0.
N_lit_fol = 0.
C_lit_frt = 0.
N_lit_frt = 0.
C_lit_crt = 0.
N_lit_crt = 0.
C_lit_tb = 0.
N_lit_tb = 0.
C_lit_stem = 0.
N_lit_stem = 0.
select case (flag_decomp)
case (20,21)
if (time .gt. 0) call read_litter_input
case(30,31)
continue
case default
! Input of litter into primary organic matter pools
! litter: x kg/tree to g/m2 (n*x*1000g/(kPatchSize m2))
! values are aggregated already as cohort
hconvd = 1000. / kpatchsize
zeig => pt%first
do while (associated(zeig))
ns = zeig%coh%species
sliti = slit(ns)
sliti%C_opm_fol = sliti%C_opm_fol + zeig%coh%litC_fol * hconvd
sliti%N_opm_fol = sliti%N_opm_fol + zeig%coh%litN_fol * hconvd
sliti%C_opm_stem = sliti%C_opm_stem + zeig%coh%litC_stem * hconvd
sliti%N_opm_stem = sliti%N_opm_stem + zeig%coh%litN_stem * hconvd
sliti%C_opm_tb = sliti%C_opm_tb + zeig%coh%litC_tb * hconvd
sliti%N_opm_tb = sliti%N_opm_tb + zeig%coh%litN_tb * hconvd
hfc = zeig%coh%litC_frt * hconvd
hfn = zeig%coh%litN_frt * hconvd
hfrtc = hconvd * zeig%coh%litC_crt
hfrtn = hconvd * zeig%coh%litN_crt
do i = 1,nroot_max
hfc1 = zeig%coh%frtrel(i)
sliti%C_opm_frt(i) = sliti%C_opm_frt(i) + hfc * hfc1
sliti%N_opm_frt(i) = sliti%N_opm_frt(i) + hfn * hfc1
sliti%C_opm_crt(i) = sliti%C_opm_crt(i) + hfrtc * hfc1
sliti%N_opm_crt(i) = sliti%N_opm_crt(i) + hfrtn * hfc1
enddo ! i (nroot_max)
C_lit_frt = C_lit_frt + zeig%coh%litC_frt
N_lit_frt = N_lit_frt + zeig%coh%litN_frt
C_lit_crt = C_lit_crt + zeig%coh%litC_crt
N_lit_crt = N_lit_crt + zeig%coh%litN_crt
C_lit_fol = C_lit_fol + zeig%coh%litC_fol
N_lit_fol = N_lit_fol + zeig%coh%litN_fol
C_lit_tb = C_lit_tb + zeig%coh%litC_tb
N_lit_tb = N_lit_tb + zeig%coh%litN_tb
C_lit_stem = C_lit_stem + zeig%coh%litC_stem
N_lit_stem = N_lit_stem + zeig%coh%litN_stem
slit(ns) = sliti
zeig => zeig%next
enddo ! show (cohorts)
do i = 1,nspec_tree
! input of delayed litter fall from dead stems
slit(i)%C_opm_tb = slit(i)%C_opm_tb + dead_wood(i)%C_tb(1)
slit(i)%N_opm_tb = slit(i)%N_opm_tb + dead_wood(i)%N_tb(1)
slit(i)%C_opm_stem = slit(i)%C_opm_stem + dead_wood(i)%C_stem(1)
slit(i)%N_opm_stem = slit(i)%N_opm_stem + dead_wood(i)%N_stem(1)
C_lit_tb = C_lit_tb + dead_wood(i)%C_tb(1)
N_lit_tb = N_lit_tb + dead_wood(i)%N_tb(1)
C_lit_stem = C_lit_stem + dead_wood(i)%C_stem(1)
N_lit_stem = N_lit_stem + dead_wood(i)%N_stem(1)
enddo ! i (nspec_tree)
! conversion g/m2/patch --> g/m2
C_lit_fol = C_lit_fol * hconvd
N_lit_fol = N_lit_fol * hconvd
C_lit_frt = C_lit_frt * hconvd
N_lit_frt = N_lit_frt * hconvd
C_lit_crt = C_lit_crt * hconvd
N_lit_crt = N_lit_crt * hconvd
C_lit_tb = C_lit_tb * hconvd
N_lit_tb = N_lit_tb * hconvd
C_lit_stem= C_lit_stem * hconvd
N_lit_stem= N_lit_stem * hconvd
end select ! flag_decomp
do j=1,nlay
cnv_opm(j) = f_cnv(C_opm(j), N_opm(j))
cnv_hum(j) = f_cnv(C_hum(j), N_hum(j))
enddo
Clitf = C_lit_frt + C_lit_crt
Nlitf = N_lit_frt + N_lit_crt
C_lit = C_lit_fol + C_lit_tb + Clitf
N_lit = N_lit_fol + N_lit_tb + Nlitf
C_lit_m = C_lit + C_lit_m
N_lit_m = N_lit + N_lit_m
C_opm = 0.
N_opm = 0.
C_opm_stem = 0.
do i = 1,nspecies
C_opm(1) = C_opm(1) + slit(i)%C_opm_frt(1) + slit(i)%C_opm_crt(1) &
+ slit(i)%C_opm_fol + slit(i)%C_opm_tb
N_opm(1) = N_opm(1) + slit(i)%N_opm_frt(1) + slit(i)%N_opm_crt(1) &
+ slit(i)%N_opm_fol + slit(i)%N_opm_tb
C_opm_stem = C_opm_stem + slit(i)%C_opm_stem
do j = 2,nlay
C_opm(j) = C_opm(j) + slit(i)%C_opm_frt(j) + slit(i)%C_opm_crt(j)
N_opm(j) = N_opm(j) + slit(i)%N_opm_frt(j) + slit(i)%N_opm_crt(j)
enddo
enddo
END subroutine cn_inp
!**************************************************************
SUBROUTINE read_litter_input
! Reading of litter input data
use data_soil_cn
use data_simul
integer lyear, lspec, ios
real helpC, helpN
logical :: lin = .TRUE.
type (species_litter) :: sliti
if (lin) read(unit_litter,*,iostat=ios) lyear, lspec, helpC, helpN
if (ios .lt. 0) lin = .FALSE.
do while (lyear .lt. time_cur)
if (lin) read(unit_litter,*,iostat=ios) lyear, lspec, helpC, helpN
if (ios .lt. 0) then
lin = .FALSE.
exit
endif
enddo
do while (lyear .eq. time_cur)
sliti = slit(lspec)
sliti%C_opm_fol = sliti%C_opm_fol + helpC
sliti%N_opm_fol = sliti%N_opm_fol + helpN
C_lit_fol = C_lit_fol + helpC
N_lit_fol = N_lit_fol + helpN
slit(lspec) = sliti
if (lin) read(unit_litter,*,iostat=ios) lyear, lspec, helpC, helpN
if (ios .lt. 0) then
lin = .FALSE.
exit
endif
enddo
if (lin) backspace (unit_litter)
END subroutine read_litter_input
!**************************************************************
SUBROUTINE n_upt(jlay)
! N uptake by roots
use data_climate
use data_par
use data_simul
use data_soil
use data_soil_cn
use help_soil_cn
use data_species
use data_stand
implicit none
! input:
integer jlay ! number of actual layer
integer i, ntr
!-----------------------------------------------------------
real NH4f, NO3f ! free available NH4-, NO3-N
real NH4u, NO3u, Nutot ! uptake of NH4-N, NO3-N, Nan_tot
real NH4jl, NO3jl ! NH4-, NO3-N
real watlay ! total water content of layer before uptake and perc.
real upt_w ! relative part of uptake water
real :: etau = 0.036 ! parameter from A. Friend (1997)
real :: fft ! temperature function of uptake from Thornley (1991)
real :: ft0 = 0. , &
ftmax = 30. , &
ftref = 20. ! parameter (°C) of temperature function from Thornley (1991)
real help, hNupt, hNupt1, Nutot1, h1, h2, N_ava, frtrel, hfrtrel, hxw
real, dimension(1:anz_coh) :: N_dem ! auxilary array for cohorts
real, external :: fred1
! no roots -> no N-uptake
if (root_fr(jlay) .lt. 1E-10) then
Nupt(jlay) = 0.
return
endif
! all NH4 and NO3 plant available
NH4jl = NH4(jlay)
NO3jl = NO3(jlay)
NH4f = NH4jl
NO3f = NO3jl
! relative part of uptake water
watlay = wats(jlay) + wupt_r(jlay)
upt_w = wupt_r(jlay) / watlay
! uptake of total available N
upt_w = 1.
fft = (temps(jlay)-ft0)*(2.*ftmax-ft0-temps(jlay))/((ftref-ft0)*(2.*ftmax-ft0-ftref))
if (fft .lt. 0.) then
fft = 0.
else
if (fft .gt. 1) fft = 1.
endif
NH4u = NH4f * fred1(jlay)
NO3u = NO3f * fft * fred1(jlay)
Nutot = (NH4u + NO3u)
Nutot1 = 0. ! actual N uptake per layer
! Uptake per cohort and m2
select case (flag_decomp)
case (0, 1, 10, 11, 20, 21, 30, 31)
if (wupt_r(jlay) .lt. 1E-10) then
Nupt(jlay) = 0.
return
else
! new balance
NH4jl = NH4(jlay) - NH4u
NO3jl = NO3(jlay) - NO3u
if (Nutot .ge. zero) then
i = 1
hxw = 0.
zeig => pt%first
do while (associated(zeig))
if (zeig%coh%species.ne.nspec_tree+2) then !exclude mistletoe
ntr = zeig%coh%ntreea
hNupt = Nutot * xwatupt(i,jlay) / wupt_r(jlay)
Nutot1 = Nutot1 + hNupt
N_ava = hNupt * kpatchsize
N_ava = N_ava /ntr ! g per tree
zeig%coh%Nuptc_d = zeig%coh%Nuptc_d + N_ava ! in g per tree
N_ava = N_ava * ntr ! Balance in g/m2
i = i+1
endif !exclusion of mistletoe
zeig => zeig%next
enddo
Nutot1 = Nutot
else
Nutot1 = 0.
do i = 1, anz_coh
xNupt(i,jlay) = 0.
enddo
endif
endif
case (40, 41)
! Ansatz A. Friend (1997, Gl. 13)
if (Nutot .ge. 1.e-6) then
i = 1
hxw = 0.
zeig => pt%first
do while (associated(zeig))
if (zeig%coh%species.ne.nspec_tree+2) then !exclude mistletoe
ntr = zeig%coh%ntreea
frtrel = zeig%coh%frtrelc(jlay) ! root percentage of the entire cohort
hNupt = frtrel * Nutot ! available nitrogen per tree cohort and layer g
hxw = hxw + frtrel
h2 = Nutot * kpatchsize
h1 = h2 * frtrel
N_ava = h1/ntr
hNupt1 = N_ava
h1 = zeig%coh%Ndemc_c - zeig%coh%Nuptc_c
h2 = zeig%coh%Ndemc_d - zeig%coh%Nuptc_d
help = h1 + h2 ! limited by actual and resudual cohort N-demand in g/m2
! limited by actual and residual cohort N-demand in g per tree
if (help .gt. N_ava) then
zeig%coh%Nuptc_d = zeig%coh%Nuptc_d + N_ava ! in g per tree
else
if (help .gt. 0.) then
zeig%coh%Nuptc_d = zeig%coh%Nuptc_d + help ! in g/m2 per Coh.
N_ava = help
else
N_ava = 0.
endif
endif
N_ava = N_ava * ntr ! balance in g/m2
h1 = N_ava
if (NH4jl .lt. h1+zero) then
h1 = h1 - NH4jl
NH4jl = zero
if (NO3jl .lt. h1+zero) then
h1 = h1 - NO3jl
NO3jl = zero
else
NO3jl = NO3jl - h1
endif
else
NH4jl = NH4jl - h1
h1 = 0.
endif
if ((NH4jl .lt. 0.) .or. (NO3jl .lt. 0.)) then
continue
endif
Nutot1 = Nutot1 + N_ava
xNupt(i,jlay) = N_ava
i = i+1
endif !exclusion of mistletoe
zeig => zeig%next
enddo
Nutot1 = Nutot1/kpatchsize
else
Nutot1 = 0.
do i = 1, anz_coh
xNupt(i,jlay) = 0.
enddo
endif
end select
NH4(jlay) = NH4jl
NO3(jlay) = NO3jl
Nupt(jlay) = Nutot1 ! N uptake per layer
END subroutine n_upt
!**************************************************************
source_code/version_2.3_windows/soil_tem.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* SOIL-Temperature - Programs *!
!* *!
!* Author: F. Suckow *!
!* *!
!* contains: *!
!* SOIL_TEMP main program for soil temperature *!
!* S_T_INI initialisation of soil temperature model *!
!* S_T_STRT initialisation of geometry parameter for the *!
!* numerical solution of the heat conduction equation *!
!* SURF_T calculation of the soil surface temperature *!
!* COND calculation of conductivity parameters *!
!* *!
!* 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_temp
! soil temperature model
use data_simul
use data_climate
use data_soil
use data_soil_t
use data_out
implicit none
integer i
! Surface temperature
call surf_t
if (flag_dayout .eq. 3) then
write (3334,*)
write (3334,*) iday
endif
! Calculation of thermal conductivity and capacity
do i=1,nlay
call cond(i,wats(i),dens(i),thick(i),pv_v(i),sandv(i),clayv(i),siltv(i),skelv(i),vol(i),spheat(i),t_cond(i),h_cap(i))
enddo ! i (nlay)
call cond(nlay, wats(nlay),dens(nlay),sh(nlay1),pv_v(nlay),sandv(nlay),clayv(nlay),siltv(nlay),skelv(nlay),vol(nlay),spheat(nlay),t_cond(nlay1),h_cap(nlay1))
call cond(nlay, wats(nlay),dens(nlay),sh(nlay2),pv_v(nlay),sandv(nlay),clayv(nlay),siltv(nlay),skelv(nlay),vol(nlay),spheat(nlay),t_cond(nlay2),h_cap(nlay2))
! Calculation of thermal diffusivity
t_cb(1) = t_cond(1)
do i=2,nlay2
t_cb(i) = (sh(i-1)*t_cond(i-1) + sh(i)*t_cond(i))/(sh(i)+sh(i-1))
enddo
if (flag_dayout .eq. 4) then
do i=1,nlay
write (3336,'(3I4, 5E11.4)') time,iday,i,watvol(i),dens(i),spheat(i),t_cond(i),h_cap(i)
enddo
write (3336, *)
endif
! Numerical solution of the heat conduction equation
call num_t
lfirst = .FALSE.
! Restore of temperature
do i=1,nlay
if (abs(sbt(i)) .lt. 1e-6) sbt(i)=0.
temps(i) = sbt(i)
enddo
! soil heat flux at soil surface
hflux_surf = 2. * t_cond(1) * (temps_surf - temps(1)) / thick(1)
1010 FORMAT (2I5, 20F8.1)
END subroutine soil_temp
!******************************************************************************
SUBROUTINE s_t_ini
! Initialisation of soil temperature model
use data_simul
use data_soil
use data_soil_t
implicit none
integer i
real, external:: kw
real tc_cont ! thermal conductivity of continuum
! Preparation of subroutine cond
! Parameter initialisation
water%tc = 0.005945 ! thermal conductivity of water at 20C J/cm/s/K
quarz%tc = 0.0879228 ! thermal conductivity of quarz at 20C
humus%tc = 0.00251 ! thermal conductivity of humus
clay%tc = 0.0251208 ! thermal conductivity of clay minerals
silt%tc = 0.02931 ! thermal conductivity of silt
air%tc = 0.00026 ! thermal conductivity of air
ice%tc = 0.021771 ! thermal conductivity of ice
stone%tc = 0.041868 ! thermal conductivity of stone
water%hc = 4.1868 ! heat capacity of water J/cm3/K
quarz%hc = 2.01 ! heat capacity of quarz
humus%hc = 2.512 ! heat capacity of humus
clay%hc = 2.01 ! heat capacity of clay minerals
silt%hc = 2.01 ! heat capacity of silt
air%hc = 0.0012 ! heat capacity of air
ice%hc = 1.884 ! heat capacity of ice
stone%hc = 1.8 ! heat capacity of stone
! shape factors
quarz%ga = 0.144 ! de Vries, S. 224
clay%ga = 0.144
silt%ga = 0.144
stone%ga = 0.144
humus%ga = 0.333
air%ga = 0.333
ice%ga = 0.125
! weighting factors for dry soil (continuous medium air)
tc_cont = air%tc
water%kwa = kw(water, tc_cont)
quarz%kwa = kw(quarz, tc_cont)
clay%kwa = kw(clay, tc_cont)
silt%kwa = kw(silt, tc_cont)
humus%kwa = kw(humus, tc_cont)
ice%kwa = kw(ice, tc_cont)
stone%kwa = kw(stone, tc_cont)
air%kwa = 1
! weighting factors for wet soil (continuous medium water)
tc_cont = water%tc
water%kww = 1
quarz%kww = kw(quarz, tc_cont)
clay%kww = kw(clay, tc_cont)
silt%kww = kw(silt, tc_cont)
humus%kww = kw(humus, tc_cont)
ice%kww = kw(ice, tc_cont)
stone%kww = kw(stone, tc_cont)
air%kww = kw(air, tc_cont)
if (flag_dayout .eq. 3) then
write (3335, '(A)') 'wet soil'
write (3335,'(6E11.4)') water%kww, air%kww, humus%kww, quarz%kww, clay%kww,ice%kww
write (3335, '(A)') 'dry soil'
write (3335,'(6E11.4)') water%kwa, air%kwa, humus%kwa, quarz%kwa, clay%kwa,ice%kwa
endif
! Calculation of thermal diffusivity
do i=1,nlay
call cond(i,wats(i),dens(i),thick(i),pv_v(i),sandv(i),clayv(i),siltv(i),skelv(i),vol(i),spheat(i),t_cond(i),h_cap(i))
enddo
call s_t_prof
call s_t_strt
! Calculation of thermal diffusivity (additional layers)
call cond(nlay, wats(nlay),dens(nlay),sh(nlay1),pv_v(nlay),sandv(nlay),clayv(nlay),siltv(nlay),skelv(nlay),vol(nlay),spheat(nlay),t_cond(nlay1),h_cap(nlay1))
call cond(nlay, wats(nlay),dens(nlay),sh(nlay2),pv_v(nlay),sandv(nlay),clayv(nlay),siltv(nlay),skelv(nlay),vol(nlay),spheat(nlay),t_cond(nlay2),h_cap(nlay2))
t_cb(1) = t_cond(1)
do i=2,nlay2
t_cb(i) = (sh(i-1)*t_cond(i-1) + sh(i)*t_cond(i))/(sh(i)+sh(i-1))
enddo
END subroutine s_t_ini
!******************************************************************************
SUBROUTINE s_t_strt
! Initialisation of geometry parameter for the
! numerical solution of the heat conduction equation
use data_soil
use data_soil_t
implicit none
integer i
real h_0, h_1
real :: ntau = 1. ! potential time step
lfirst = .TRUE.
nlay1 = nlay+1
nlay2 = nlay+2
sh(1) = thick(1)
sb(1) = 2. / sh(1)
sv(mfirst) = sh(mfirst)
sbt(mfirst) = temps_surf
do i=mfirst+1,nlay
sbt(i) = temps(i)
sh(i) = thick(i)
enddo
sbt(nlay1) = temps(nlay)
sbt(nlay2) = temps(nlay)
sh(nlay1) = 2. * thick(nlay)
sh(nlay2) = 100.
h_0 = sh(1)
do i= mfirst+1, nlay2
h_1 = sh(i)
sb(i) = 2. / (h_1 + h_0)
sv(i) = h_1 * ntau
h_0 = h_1
enddo
END subroutine s_t_strt
!******************************************************************************
SUBROUTINE surf_t
! Calculation of soil surface temperature
use data_climate
use data_simul
use data_soil
use data_soil_t
use data_stand
implicit none
real day
real cof ! daily correction cefficient
real dampcof ! stand damping coefficient
real helplai ! thermal conductivity of organic layer (global vereinbaren und vom Vortag merken!!!)
integer unit_tmp, helptyp
character(80) text
! read surface temperature; Oberflaechentemperatur einlesen
if (flag_surf .eq. 2) then
if (lfirst) then
write (*,'(A)', advance='no') 'Reading of soil surface temperature, please type file name:'
read (*,'(A)') text
unit_tmp = getunit()
open (unit_tmp, file=trim(text), status='unknown')
read (unit_tmp,'(A)') text
read (unit_tmp, *) day, temps_surf
return
else
read (unit_tmp, *) day, temps_surf
return
endif
endif
! snow
if (snow .lt. 0.05) then ! calculation of temps_surf in subroutine snowpack
dampcof = 1.0
if (waldtyp .ge. 110 .and. (waldtyp .ne. 125)) then
helptyp = 110
else
helptyp = waldtyp
endif
select case (helptyp)
case (10,20,25,30,31,35,37,38,70,71,75,76,125) ! Spruce; Fichte
if (iday .lt. 90 .or. (iday .gt. 320)) then
dampcof=0.8
else if (iday .lt. 115) then
dampcof=1.0
else if (iday .gt. 240) then
dampcof=1.0
else
dampcof=0.7
endif
case (40,50,51,52,54,55,56,60,61,62,64,65,66,90,100) ! Pine; Kiefer
if (iday .lt. 90 .or. (iday .gt. 320)) then
dampcof=1.5
else if (iday .lt. 115) then
dampcof=1.2
else if (iday .gt. 285) then
dampcof=1.3
else
dampcof=0.8
endif
case (110) ! Beech and other decidous trees; Buche und andere Laubhoelzer
if (LAI .gt. 1.) then
if (iday .gt. 50) then
if (iday .lt. 100 .or. (iday .gt. 300 .and. iday .lt. 345)) then
dampcof=1.2
else if (iday .gt. 130 .and. iday .le. 300) then ! for beech; fuer Buche
dampcof=0.8 ! for beech; fuer Buche
endif
endif
else
dampcof=1.2 ! for beech; fuer Buche
endif
end select
! Daempfung berechnen nach Paul et al. (2004)
day = iday
cof = abs(-0.00003*day*day + 0.0118*day - 0.0703)
if (flag_surf .eq. 0) then
temps_surf = (c0*airtemp + c1*airtemp_1 + c2*airtemp_2) * cof * dampcof
temps(1) = temps_surf
else
if (flag_surf .eq. 3) then
cof = 1
dampcof = 1.0
endif
temps_surf = (c0*airtemp + c1*airtemp_1 + c2*airtemp_2) * cof * dampcof
endif
endif ! snow
if (flag_dayout .eq. 3) then
write (1222,'(A,I5,F10.4,3F8.2)') 'day, cof, dampcof', iday, cof, dampcof, temps_surf, airtemp
endif
END subroutine surf_t
!******************************************************************************
SUBROUTINE cond(ilay,watsi,densi,thicki,pvi,sandi,clayi,silti,skelvi,voli,spheati,tcondi,hcapi)
! Calculation of thermal conductivity and capacity
! de Vries-approach
use data_par
use data_soil
use data_soil_cn
use data_soil_t
use data_simul
implicit none
! input
integer ilay ! number of layer
real watsi ! water content mm
real densi ! soil density
real thicki ! layer thickness
real spheati ! specific heat capacity
real dmi ! dry mass g/m2
real pvi ! pore volume
real quarzi ! quarz fraction in soild soil
real sandi ! sand fraction in soild soil
real clayi ! clay fraction in soild soil
real silti ! silt fraction in soild soil
real skelvi ! skeleton fraction in soil
real tc_cont ! thermal conductivity of continuum
real wcvol ! water content (vol%)
! output
real tcondi, tcond0, tcond1, tcond2, tcond3 ! thermal conductivity
real hcapi, hcap0, hcap1, hcap2, hcap3 ! thermal capacity
real numera, denom ! numerator, denominator of calculation of thermal conductivity
real hum_dens, densi1, pvi1, hvf, hvf1
real aa, bb, cc, dd, vfm, vfs, massfr ! Campbell-Ansatz
real wkw,akw,hkw,qkw,ckw,skw,ikw,tkw, skel, voli, restvol
! density g/cm3
hum_dens = 1.3 !Density of humus (compressed, without air)
quarzi = sandi
! dry mass
dmi = voli * densi
voli = thicki * 10000.
hvf = (C_opm(ilay) + C_hum(ilay)) / cpart ! Masse (g)
! volume fractions
skel = 1. - skelvi
pvi1 = skel * pvi/100.
water%vf = skel * watsi/(10.*thicki)
air%vf = pvi1 - water%vf
if (air%vf .lt. 0.) then
continue
endif
hvf = hvf / hum_dens ! volume; Volumen
restvol = voli - (skelvi + pvi1)*voli - hvf
humus%vf = hvf / voli
quarz%vf = quarzi*restvol / voli
clay%vf = clayi*restvol / voli
silt%vf = silti*restvol / voli
stone%vf = skelvi
ice%vf = 0.
if (flag_dayout .ge. 3) then
write (3334,'(3I4,F8.3,8F10.4)') time,iday,ilay,pvi1, water%vf, air%vf,humus%vf,quarz%vf,clay%vf,silt%vf,stone%vf,ice%vf
if (ilay .eq. nlay) write (3334, *)
endif
select CASE (flag_cond)
CASE (1, 11, 21, 31, 41) ! Neusypina
if (densi .lt. 0.6) then
densi1 = 0.6
else
densi1 = densi
endif
wcvol = watsi/(10.*thicki)
tcondi = ((3.*densi1-1.7)*0.001)/(1.+(11.5-5.*densi1) &
*EXP((-50.)*(wcvol/densi1)**1.5))*86400.
tcondi = tcondi * 4.1868 ! convertation cal/(cm s K) in J/(cm s K)
! heat capacity J/(cm3 K)
hcapi = densi1*spheati + wcvol*4.1868
hcap1 = hcapi
tcond1 = tcondi
CASE (0, 10, 20, 30, 40) ! de Vries
! Determination of continuous medium
if (watsi .gt. 0.95 * pv(ilay)) then
! wet soil
wkw = water%kww
akw = air%kww
hkw = humus%kww
qkw = quarz%kww
ckw = clay%kww
skw = silt%kww
tkw = stone%kww
ikw = ice%kww
else
! dry soil
wkw = water%kwa
akw = air%kwa
hkw = humus%kwa
qkw = quarz%kwa
ckw = clay%kwa
skw = silt%kwa
tkw = stone%kwa
ikw = ice%kwa
endif
numera = wkw * water%vf * water%tc + qkw * quarz%vf * quarz%tc + ckw * clay%vf * clay%tc + &
skw * silt%vf * silt%tc + hkw * humus%vf * humus%tc + akw * air%vf * air%tc + &
tkw * stone%vf * stone%tc + ikw * ice%vf * ice%tc
denom = wkw * water%vf + qkw * quarz%vf + ckw * clay%vf + skw * silt%vf + &
hkw * humus%vf + akw * air%vf + tkw * stone%vf + ikw * ice%vf
tcond0 = numera/denom * 86400. ! s --> day
CASE(3, 13, 23, 33, 43) ! Campbell
vfm = clay%vf + silt%vf + stone%vf
vfs = vfm + quarz%vf + humus%vf
if (watsi .gt. 0.95 * pv(ilay)) then
! wet soil
aa = 0.57 + 1.73*quarz%vf + 0.93*vfm
aa = aa / (1. - 0.74*quarz%vf - 0.49*vfm) - 2.8*vfs*(1.-vfs)
bb = 2.8 * vfs * water%vf
tcond3 = (aa + bb * water%vf) ! W/m/K
else if (watsi .le. wilt_p(ilay)) then
! dry soil
tcond3 = 0.03 + 0.7 * vfs * vfs ! W/m/K
else
massfr = 2.65 * (vfm + quarz%vf) + 1.3 * humus%vf
massfr = 2.65 * clay%vf / massfr
aa = 0.57 + 1.73*quarz%vf + 0.93*vfm
aa = aa / (1. - 0.74*quarz%vf - 0.49*vfm) - 2.8*vfs*(1.-vfs)
bb = 2.8 * vfs * water%vf
cc = 1. + 2.6 * sqrt(clay%vf)
dd = 0.03 + 0.7 * vfs * vfs
tcond3 = aa + bb*water%vf - (aa-dd) * exp(-(cc*water%vf)**4) ! W/m/K
endif
tcond3 = tcond3 / 100. ! W/m/K ==> J/(cm s K)
tcond3 = tcond3 * 86400. ! s --> day
end select
! heat capacity J/(cm3 K)
hcap0 = water%vf * water%hc + quarz%vf * quarz%hc + clay%vf * clay%hc + silt%vf * silt%hc + &
humus%vf * humus%hc + air%vf * air%hc + stone%vf * stone%hc + ice%vf * ice%hc
if (flag_dayout .eq. 4) then
write (3337,'(3I4, 6E11.4)') time,iday,ilay,tcond0,tcond1,tcond2,tcond3,hcap0,hcap1
if (ilay .eq. nlay) write (3337, *)
endif
select CASE (flag_cond)
CASE (0, 10, 20, 30, 40) ! de Vries
hcapi = hcap0
tcondi = tcond0
CASE (1, 11, 21, 31, 41) ! Neusypina
hcapi = hcap1
tcondi = tcond1
CASE (3, 13, 23, 33, 43) ! Campbell
hcapi = hcap0
tcondi = tcond3
end select
END subroutine cond
!**************************************************************
real FUNCTION kw(part, tc_cont)
! Function for calculating weighting factor k
! in calculating thermal conductivity
use data_soil_t
implicit none
type (therm_par):: part ! soil fraction (particles)
real tc_cont ! thermal conductivity of continuum
real term, ga
term = part%tc / tc_cont -1.
ga = part%ga
kw = (2./(1.+ term*ga) + 1./(1.+ term*(1.-2.*ga)))/3.
end FUNCTION
!******************************************************************************
SUBROUTINE num_t
! Numerical solution of the heat conduction equation
use data_soil
use data_soil_t
use data_simul
implicit none
integer i
logical lcase ! logical control of Cholesky procedure
real hflux ! heat flux at surface (right side)
lcase = .TRUE.
! Determination of the volume matrix
svv = sv * h_cap
if (lfirst) svva = svv
! Determination (side diagonal; Nebendiagonale) !
do i=1,nlay2
son(i) = -sb(i) * t_cb(i)
enddo
son(nlay2+1) = 0.0
! Determination (main diagonal; Hauptdiagonale) !
do i=1,nlay2
soh(i) = svv(i) - son(i) - son(i+1)
enddo
hflux = temps_surf * sb(1) * t_cb(1) ! Set heat flux at surface at right side
if (.not.lfirst) then
! Calculation of the right side
do i=1,nlay2
sxx(i) = (svva(i) + (svv(i)-svva(i))/sh(i)) * sbt(i)
enddo
sxx(1) = sxx(1) + hflux
! Iteration (Cholesky procedure)
call chl3 (nlay2, son, soh, sxx, lcase)
! Results of iteration on temperature help array
sbt = sxx
endif ! lfirst
! Restore of geometry matrix
svva = svv
END subroutine num_t
!******************************************************************************
SUBROUTINE chl3 (n, a, b, x, lcase)
! Solution of EX = Z (E - tridiagonal, symmetric matrix)
! with Cholesky procedure (E = LDL')
implicit none
! input
integer n ! rang of matrix
logical lcase ! logical control of Cholesky procedure
! .TRUE. for start of iteration
real, dimension(n) :: a, & ! Nebendiagonale
b ! main diagonal
! output
real, dimension(n) :: x ! solution vector
! local variables
integer i, j, j1
real, dimension(n) :: d, ul
! Calculation of the left upper triangle matrix L
! and of the diagonal matrix D at start of iteration
if (lcase) then
d(1) = b(1)
do i=2,n
ul(i) = a(i) / d(i-1)
d(i) = b(i) - ul(i)*a(i)
enddo
lcase = .FALSE.
endif
! Solution of LY = Z
do i=2,n
x(i) = x(i) - ul(i)*x(i-1)
enddo
! Solution of L'X = D(-1)Y
x(n) = x(n) / d(n)
do i=1,n-1
j = n-i
j1 = j+1
x(j) = x(j)/d(j) - ul(j1)*x(j1)
enddo
END subroutine chl3
!******************************************************************************
source_code/version_2.3_windows/soil_tem_ini.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* 4C (FORSEE) Simulation Model *!
!* *!
!* *!
!* contains: *!
!* s_t_prof generates initial soil temp. profile *!
!* BTFOUR TRICOF use to develope soil-surface-temp. *!
!* *!
!* 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 s_t_prof
! Generation of initial soil temperature profile
use data_par
use data_soil
use data_soil_t
use data_simul
implicit none
integer i, ia, ie, k
real ath, dth, rfr, sn, uhf, u, vk, vh, fourterm
real tcsu, hcsu
ia = 1
ie = 365
TQ = 10.
CALL BTFOUR
tcsu = 0.
hcsu = 0.
rfr = 2. * pi / 365. ! radial frequency; Radialfrequenz
UHF=2.*pi/(IE-IA+1)
u = uhf * it
! calculation of temperature profile commonly from day 1 (it=1) set in data_soil;
! Temperaturprofil berechnen, standardmaessig it=1 (1. Tag) in data_soil gesetzt
do i = 1, nlay
tcsu = tcsu + t_cond(i)*thick(i)
hcsu = hcsu + h_cap(i)*thick(i)
ath = tcsu / hcsu ! for a weighted mean both values are divided by the depth (i), thus they cancel each other; fuer gewichtetes Mittel beide Werte durch depth(i) teilen ==> weggekuerzt
DTH=SQRT(2*ATH/RFR)
VH=mid(I)/DTH
fourterm = 0.
do k = 1, nk
VK=VH*SQRT(K+0.)
SN=FTA(K)*EXP(-VK)*SIN(U*K+FTO(K)-VK)
fourterm = fourterm + SN
enddo
temps(i) = TQ + fourterm
if (flag_dayout .eq. 3) write (2244, *) i, temps(i), mid(i), ath, dth, fourterm
enddo
END subroutine s_t_prof
!******************************************************************************
SUBROUTINE BTFOUR
! using TRICOF for a Fourier series development for ground surface temperature;
! Fourierreihenentwicklung fuer Boden-Oberflaechen-Temperatur unter Nutzung von TRICOF
use data_climate
use data_par
use data_soil
use data_soil_t
use data_simul
implicit none
integer i, n, nt, nts, nf, nend, naf, no, ne, lf
real a0
real, dimension(184):: FA,FB
! set amount of auxiliary points NF for transformation;
! Anzahl der Stuetzstellen NF fuer Transformation festlegen
nend = 365
naf = 1
NT=NEND-NAF+1
NTS=1
NF=(NT+NTS-1)/NTS
N=(NF-1)/2
IF((2*N-NF+1) .LT. 0) THEN
NF=NF-1
NT=(NF*NTS)-NTS+1
NEND=NAF+NT-1
ENDIF
NE=1+(NF/2)
NO=NE-2
NK=NO
! calculation of auxiliary points; Stuetzstellen berechnen
LF=NAF
DO I=1,NF
airtemp = tp(lf,1)
airtemp_1 = tp(lf-1,1)
airtemp_2 = tp(lf-2,1)
rad = rd(lf,1)
iday = lf
call surf_t
if (lf .eq. 1) temps(1) = temps_surf
Four_sp(i) = temps_surf
LF=LF+NTS
ENDDO
! Fourier transformation; FOURIERTRANSFORMATION
CALL TRICOF(Four_sp,NF,FA,NE,FB,NO,1)
A0 = FA(1) / 2.
TQ = A0
! coefficient to transform solution; Koeffizienten fuer Loesung transformieren
DO I=1,NK
FTA(I) = SQRT(FA(I+1)*FA(I+1) + FB(I)*FB(I))
FTA(I) = FTA(I) * SIGN(1.,FB(I))
if(FB(I).eq. 0.) then
FTO(I) = pi/2.
else
FTO(I) = ATAN(FA(I+1)/FB(I))
end if
FTO(I) = FTO(I) - (NEND+NAF)*PI*I/(NEND-NAF)
ENDDO
END SUBROUTINE BTFOUR
!******************************************************************************
source_code/version_2.3_windows/sorting.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* Subroutines for: *!
!* - dimsort: sorting of cohorts according to a *!
!* charcteristic variable *!
!* - sort2: subroutine from Numerical recipes *!
!* *!
!* 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 dimsort(n,var,ranktable)
USE data_species
USE data_stand ! state variables of stand, cohort and cohort element
IMPLICIT NONE
INTEGER :: isort,n
INTEGER :: ranktable(n)
CHARACTER(3) :: var
REAL :: sortarray(n)
ranktable=0
sortarray=0
isort=1
zeig=>pt%first
DO
IF (.not.ASSOCIATED(zeig)) exit
IF (zeig%coh%species .le. nspec_tree) THEN ! for trees only
ranktable(isort) = zeig%coh%ident
IF(var=='hei') sortarray(isort) = zeig%coh%height
IF(var=='dbh') sortarray(isort) = zeig%coh%diam
isort=isort+1
ENDIF
zeig=>zeig%next
END DO
CALL sort2(n,sortarray,ranktable)
END SUBROUTINE dimsort
!******************************************************************************
SUBROUTINE sort2(n,arr,brr)
! sorts array arr(1:n) into an ascending order and
! makes the corresponding rearrangement of the array brr(1:n)
INTEGER n,M,NSTACK
REAL arr(n)
INTEGER brr(n)
PARAMETER (M=7,NSTACK=50)
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
REAL a,b,temp
jstack=0
l=1
ir=n
1 if(ir-l.lt.M)then
do 12 j=l+1,ir
a=arr(j)
b=brr(j)
do 11 i=j-1,1,-1
if(arr(i).le.a)goto 2
arr(i+1)=arr(i)
brr(i+1)=brr(i)
11 continue
i=0
2 arr(i+1)=a
brr(i+1)=b
12 continue
if(jstack.eq.0)return
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
k=(l+ir)/2
temp=arr(k)
arr(k)=arr(l+1)
arr(l+1)=temp
temp=brr(k)
brr(k)=brr(l+1)
brr(l+1)=temp
if(arr(l+1).gt.arr(ir))then
temp=arr(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
temp=brr(l+1)
brr(l+1)=brr(ir)
brr(ir)=temp
endif
if(arr(l).gt.arr(ir))then
temp=arr(l)
arr(l)=arr(ir)
arr(ir)=temp
temp=brr(l)
brr(l)=brr(ir)
brr(ir)=temp
endif
if(arr(l+1).gt.arr(l))then
temp=arr(l+1)
arr(l+1)=arr(l)
arr(l)=temp
temp=brr(l+1)
brr(l+1)=brr(l)
brr(l)=temp
endif
i=l+1
j=ir
a=arr(l)
b=brr(l)
3 continue
i=i+1
if(arr(i).lt.a)goto 3
4 continue
j=j-1
if(arr(j).gt.a)goto 4
if(j.lt.i)goto 5
temp=arr(i)
arr(i)=arr(j)
arr(j)=temp
temp=brr(i)
brr(i)=brr(j)
brr(j)=temp
goto 3
5 arr(l)=arr(j)
arr(j)=a
brr(l)=brr(j)
brr(j)=b
jstack=jstack+2
if(jstack.gt.NSTACK)pause 'NSTACK too small in sort2'
if(ir-i+1.ge.j-l)then
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif
endif
goto 1
END Subroutine
! (C) Copr. 1986-92 Numerical Recipes Software "!D#+.
source_code/version_2.3_windows/stand_bal.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* Subroutines for: *!
!* - STAND_BALANCE: Recalculation of stand variables *!
!* contains: *!
!* UPDATE_AGE *!
!* - STAND_BAL_SPEC *!
!* - CLASS *!
!* - CLASST *!
!* - CLASS_MAN *!
!* - CALC_HEIDOM *!
!* - MAX_HEIGHT(nrmax,anz,cohl) *!
!* - STANDUP: Update of cover and ceppot *!
!* - LITTER: Summation variables of litter fractions *!
!* - calc_ind_rep: calculation of representation index *!
!* - overstorey *!
!* *!
!* 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 stand_balance
use data_species
use data_stand
use data_climate
use data_simul
use data_manag
use data_out
use data_par
implicit none
integer i, ntr, nd, hanz
integer, dimension(nspecies) :: helpin, helpout
real conv ! conversion factor
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time, ' stand_balance'
if(time>0. .and. flag_standup.ne.2) call update_age
! calculation of total dead biomass per cohort and total biomass of allcohorts
! calc. of ceppot
anz_sveg = 0
anz_tree = 0.
anz_tree_in = 0.
anz_tree_out = 0.
anz_spec_in = 0.
anz_spec_out = 0.
anz_coh_in = 0.
anz_coh_out = 0.
anz_coh_act = 0.
lai_in = 0.
lai_out = 0.
totfol_in = 0.
totfol_out = 0.
med_diam_in = 0.
med_diam_out = 0.
hmean_in = 0.
hmean_out = 0.
mean_height = 0.
sumbio = 0.
sumbio_in = 0.
sumbio_out = 0.
sumNPP = 0.
drIndAl = 0.
Ndem = 0.
helpin = 0
helpout = 0
basal_area = 0.
totstem_m3 = 0.
totsteminc_m3 = 0.
totsteminc = 0
autresp = 0.
totfol = 0.
totsap = 0.
totfrt = 0.
totfrt_p = 0.
totcrt = 0.
tottb = 0.
tothrt = 0.
sumbio_sv = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
ns = zeig%coh%species
ntr = zeig%coh%ntreeA
svar(ns)%daybb = zeig%coh%day_bb
if(ns.le.nspec_tree) then
if(zeig%coh%ident .le. coh_ident_max) then
anz_coh_act = anz_coh_act + 1
anz_tree = anz_tree + ntr
zeig%coh%totBio = zeig%coh%x_fol + zeig%coh%x_sap + zeig%coh%x_hrt + zeig%coh%x_tb + zeig%coh%x_frt +zeig%coh%x_crt
zeig%coh%Dbio = zeig%coh%nTreeD * zeig%coh%totBio
sumbio = sumbio + ntr * zeig%coh%totBio
sumNPP = sumNPP + ntr * zeig%coh%NPP
Ndem = Ndem + ntr * zeig%coh%Ndemc_c
select case (flag_dis)
case (0,1)
autresp = autresp + ntr * zeig%coh%maintres
case (2)
autresp = autresp + ntr * (zeig%coh%maintres+zeig%coh%biocost_all)
end select
totfol = totfol + ntr * zeig%coh%x_fol
totsap = totsap + ntr * zeig%coh%x_sap
totfrt = totfrt + ntr * zeig%coh%x_frt
totcrt = totcrt + ntr * zeig%coh%x_crt
tottb = tottb + ntr * zeig%coh%x_tb
tothrt = tothrt + ntr * zeig%coh%x_hrt
if (zeig%coh%height.le.thr_height) then
seedlfrt = seedlfrt + zeig%coh%x_frt * ntr
endif
totstem_m3 = totstem_m3 + (ntr*zeig%coh%x_sap + ntr*zeig%coh%x_hrt) &
/(spar(ns)%prhos*1000000) ! conversion kg/patch ---m³/ha
nd = zeig%coh%nDaysGr
if (nd .gt. 0) drIndAl = drIndAl + ntr * zeig%coh%drIndAl * zeig%coh%NPP / nd
endif
if(zeig%coh%ident > coh_ident_max) then
anz_tree_in = anz_tree_in + ntr
sumbio_in = sumbio_in + ntr * zeig%coh%totBio
anz_coh_in = anz_coh_in + 1
helpin(ns) = ns
lai_in = lai_in + ntr * zeig%coh%t_leaf/kpatchsize
totfol_in = totfol_in + ntr * zeig%coh%x_fol
med_diam_in = med_diam_in + ntr * (zeig%coh%diam**2)
hmean_in = hmean_in + ntr * zeig%coh%height
totfrt = totfrt + ntr * zeig%coh%x_frt
endif
if((zeig%coh%nTreeD > 0.1) .or. (zeig%coh%nTreeM > 0.1) .or. (zeig%coh%nTreet > 0.1)) then
hanz = zeig%coh%nTreeD + zeig%coh%nTreeM + zeig%coh%nTreet
anz_tree_out = anz_tree_out + hanz
sumbio_out = sumbio_out + hanz * zeig%coh%totBio
sumNPP = sumNPP + hanz * zeig%coh%NPP ! eliminated (died or harvested) trees produce during the year as well;
autresp = autresp + hanz * zeig%coh%maintres
anz_coh_out = anz_coh_out + 1
helpout(ns) = ns
lai_out = lai_out + hanz * zeig%coh%t_leaf/kpatchsize
totfol_out = totfol_out + hanz * zeig%coh%x_fol
med_diam_out = med_diam_out + hanz * (zeig%coh%diam**2)
hmean_out = hmean_out + hanz * zeig%coh%height
endif
else
ntr = zeig%coh%ntreeA
anz_sveg = anz_sveg +1
zeig%coh%totBio = zeig%coh%x_fol + (1.+spar(ns)%alphac)*(zeig%coh%x_sap + zeig%coh%x_hrt) + zeig%coh%x_frt
sumbio_sv = sumbio_sv + ntr * zeig%coh%totBio
totfrt_p = totfrt_p + ntr * zeig%coh%x_frt
end if !only trees
zeig=>zeig%next
end do
if (flag_cumNPP .eq. 1) then
cum_sumNPP = cum_sumNPP + sumNPP
flag_cumNPP = 0
endif
if (sumNPP .gt. 1E-06) drIndAl = drIndAl / sumNPP
! conversion kg/patch ---> kg/ha; N/patch ---> N/ha
conv = 10000./kpatchsize
totfrt_p = totfrt_p + totfrt ! Rootmass f. patch (trees and soil veg.) save before conversion; Wurzelmenge vor Umrechnung sichern
if (totfrt_p .gt. zero) then
totfrt_1 = 1./totfrt_p ! reciprocal for later calculationshKehrwert f. spaetere Berechnungen
else
totfrt_1 = 0.
endif
totfrt = totfrt * conv
totfol = totfol * conv
totfol_in = totfol_in * conv
totfol_out = totfol_out * conv
tottb = tottb * conv
totsap = totsap * conv
tothrt = tothrt * conv
totcrt = totcrt * conv
sumbio = sumbio * conv
sumbio_in = sumbio_in * conv
sumbio_out = sumbio_out * conv
sumbio_sv = sumbio_sv * conv
Ndem = Ndem / kpatchsize ! g per tree --> g/m2
totstem_m3 = totstem_m3* conv ! m3/ha
anz_tree_ha = anz_tree * conv
anz_tree_in = anz_tree_in * conv
anz_tree_out = anz_tree_out * conv
do i=1, nspec_tree+1 ! for all but mistletoe
if (helpin(i) > 0) anz_spec_in = anz_spec_in + 1
if (helpout(i) > 0) anz_spec_out = anz_spec_out + 1
enddo
if(anz_tree_in > 0.) then
med_diam_in = sqrt(med_diam_in/anz_tree_in)
hmean_in = hmean_in / anz_tree_in
endif
if(anz_tree_out > 0.) then
med_diam_out = sqrt(med_diam_out/anz_tree_out)
hmean_out = hmean_out / anz_tree_out
endif
! call species values
call classt
call stand_bal_spec
call calc_ind_rep
!call classification of stand diameter and height
call class
! moving of understorey tree cohorts to overstorey tree cohorts
if(flag_mg.ne.33) call overstorey
contains
!**************************************************************
subroutine update_age
if(flag_standup.ne. 2) then
zeig=>pt%first
do
if(.not.associated(zeig)) exit
zeig%coh%x_age=zeig%coh%x_age + 1
zeig=>zeig%next
end do
end if
end subroutine update_age
end subroutine stand_balance
!**************************************************************
subroutine stand_bal_spec
use data_climate
use data_out
use data_simul
use data_site
use data_stand
use data_species
use data_par
use data_manag
implicit none
integer :: i, j, k, ntr, nd, lowtree, hntr, spec_new
real,dimension(nspec_tree):: vgldom1, vgldom2, vgldom_spec1, vgldom_spec2
integer,dimension(nspec_tree):: anzdom1, anzdom2, anzdom_spec1, anzdom_spec2, &
helpdiam
integer,dimension(nspecies):: helpanz
integer helpntr
integer help_nr_inf_trees
logical lhelp
INTEGER leapyear
real atemp, hh, help_height_top
real triangle
real, external :: daylength
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time, ' stand_bal_spec'
! Initialisation
vgldom1 = 0.
vgldom2 = 0.
anzdom1 = 0
anzdom2 = 0
med_diam = 0.
mean_diam = 0.
mean_height = 0.
hdom = 0. ! dominante height (highest or two highest cohorst); Hoehe (hoehste oder die beiden hoechsten Kohorten)
anz_spec = 0 ! currently existing species
lowtree = 0 ! amount of trees with DBH=0 for the whole population; Anzahl Baeume mit DBH = 0 fuer gesamten Bestand
hntr = 0
helpanz = 0 ! auxiliary variable to count species; Hilfsvariable um Spezies zu zaehlen
helpdiam = 0 ! amount of trees with DBH=0 per species; Anzahl Baeume mit DBH = 0 pro Spezies
vgldom_spec1 = 0.
vgldom_spec2 = 0.
anzdom_spec1 = 0
anzdom_spec2 = 0
svar%med_diam = 0.
svar%mean_diam = 0.
svar%mean_jrb = 0.
svar%mean_height= 0.
svar%basal_area = 0.
svar%sumNPP = 0.
svar%Ndem = 0.
svar%Nupt = 0.
svar%sum_ntreea = 0.
svar%sum_ntreed = 0.
svar%sum_bio = 0.
svar%sum_lai = 0.
svar%anz_coh = 0.
svar%drIndAl = 0.
svar%totsteminc = 0.
svar%totsteminc_m3 = 0.
svar%fol = 0.
svar%sap = 0.
svar%hrt = 0.
svar%frt = 0.
!update height of mistletoe to height of mistletoe-infected cohort
if (flag_mistle.ne.0) then
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
if (zeig%coh%mistletoe.eq.1) then
help_height_top=zeig%coh%height
end if
zeig => zeig%next
ENDDO
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
if (zeig%coh%species.eq.nspec_tree+2) then
zeig%coh%height = help_height_top !upper crown
zeig%coh%x_hbole = zeig%coh%height-50. !lower crown
end if
zeig => zeig%next
ENDDO
end if ! end update of height of Mistletoe
! update of #of mistletoe upon dist_manag
if (flag_mistle.ne.0) then
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
if (zeig%coh%mistletoe.eq.1) then
help_nr_inf_trees=zeig%coh%nTreeA
end if
zeig => zeig%next
ENDDO
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
if (zeig%coh%species.eq.nspec_tree+2) then
zeig%coh%nTreeA= help_nr_inf_trees*AMAX1(1.,dis_rel(time))
zeig%coh%nta=zeig%coh%nTreeA
end if
zeig => zeig%next
ENDDO
end if ! end update #of mistletoe
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
ns = zeig%coh%species
helpanz(ns) = helpanz(ns) + 1 ! all species incl. ground vegetation;
IF(zeig%coh%ident .le. coh_ident_max) THEN
ntr = zeig%coh%ntreea
IF(ns .le. nspec_tree) THEN
!
IF((ns .le. nspec_tree) .and. (ntr > 0.) .and. (zeig%coh%diam > 0.)) THEN
svar(ns)%med_diam = svar(ns)%med_diam + ntr * (zeig%coh%diam**2)
med_diam = med_diam + ntr * (zeig%coh%diam**2)
mean_diam = mean_diam + ntr*zeig%coh%diam
svar(ns)%mean_diam = svar(ns)%mean_diam + ntr*zeig%coh%diam
svar(ns)%mean_height = svar(ns)%mean_height + ntr*zeig%coh%height
svar(ns)%mean_jrb = svar(ns)%mean_jrb + ntr*zeig%coh%jrb
mean_height = mean_height + ntr*zeig%coh%height
hntr = hntr + ntr
ELSE
! Trees with DBH=0 for population and per species; Baeume mit DBH =0 fuer Bestand und pro Spezies
lowtree = lowtree + ntr
helpdiam(ns) = helpdiam(ns) + ntr
ENDIF ! ns
IF(zeig%coh%height > vgldom1(ns)) THEN
vgldom2(ns) = vgldom1(ns)
anzdom2(ns) = anzdom1(ns)
vgldom1(ns) = zeig%coh%height
anzdom1(ns) = ntr
ELSE
if(zeig%coh%height > vgldom2(ns))then
vgldom2(ns) = zeig%coh%height
anzdom2(ns) = ntr
endif
ENDIF ! vgldom1
IF(zeig%coh%height > vgldom_spec1(ns)) THEN
vgldom_spec2(ns) = vgldom_spec1(ns)
anzdom_spec2(ns) = anzdom_spec1(ns)
vgldom_spec1(ns) = zeig%coh%height
anzdom_spec1(ns) = ntr
ELSE
if(zeig%coh%height > vgldom_spec2(ns))then
vgldom_spec2(ns) = zeig%coh%height
anzdom_spec2(ns) = ntr
endif
ENDIF ! vgldom_spec2
ELSE
svar(ns)%dom_height = zeig%coh%height
ENDIF ! end loop across trees
svar(ns)%sumNPP = svar(ns)%sumNPP + ntr * zeig%coh%NPP
svar(ns)%sum_ntreea = svar(ns)%sum_ntreea + ntr
svar(ns)%sum_ntreed = svar(ns)%sum_ntreed + zeig%coh%nTreeD + zeig%coh%nTreeM ! died or harvested trees of current year; ausgeschiedene Bäume des akt. Jahres
svar(ns)%Ndem = svar(ns)%Ndem + ntr * zeig%coh%Ndemc_c
svar(ns)%Nupt = svar(ns)%Nupt + ntr * zeig%coh%Nuptc_c
svar(ns)%sum_bio = svar(ns)%sum_bio + ntr * zeig%coh%totBio
svar(ns)%sum_lai = svar(ns)%sum_lai + ntr * zeig%coh%t_leaf/kpatchsize
svar(ns)%anz_coh = svar(ns)%anz_coh + 1
svar(ns)%totsteminc = svar(ns)%totsteminc + ntr * zeig%coh%stem_inc
if (zeig%coh%species.ne.nspec_tree+2) then !no stem increment for mistletoe
svar(ns)%totsteminc_m3 = svar(ns)%totsteminc_m3 + ntr * zeig%coh%stem_inc /(spar(ns)%prhos*1000000)
endif
svar(ns)%fol = svar(ns)%fol + ntr * zeig%coh%x_fol
svar(ns)%sap = svar(ns)%sap + ntr * zeig%coh%x_sap
svar(ns)%hrt = svar(ns)%hrt + ntr * zeig%coh%x_hrt
svar(ns)%frt = svar(ns)%frt + ntr * zeig%coh%x_frt
nd = zeig%coh%nDaysGr
if (nd .gt. 0) svar(ns)%drIndAl = svar(ns)%drIndAl + ntr * zeig%coh%drIndAl * zeig%coh%NPP / nd
ENDIF ! coh%ident
zeig%coh%ntreed = 0.
zeig%coh%ntreem = 0.
zeig => zeig%next
ENDDO ! cohort loop
! neue Spezies feststellen und belegen
if (time .gt. 1) then
do i=1,nspecies
if (helpanz(i) > 0) then
spec_new = 0
lhelp = .True.
do j=1,anrspec
if (nrspec(j) .eq. i) lhelp = .False.
enddo
if (lhelp) then
spec_new = i
if(spec_new.le.nspec_tree) then
IF(spar(spec_new)%Phmodel==1) THEN
svar(spec_new)%Pro = 0.
svar(spec_new)%Inh = 1.
ELSE
svar(spec_new)%Pro = 0.
svar(spec_new)%Inh = 0.
svar(spec_new)%Tcrit = 0.
END IF
! initialize pheno state variables with climate from the actual year
do j = spar(ns)%end_bb+1, 365
atemp = tp(j, time)
hh = DAYLENGTH(j,lat)
SELECT CASE(ns)
CASE(1,8)
!Fagus
! Promotor-Inhibitor model 11
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa* &
triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,atemp)* &
(1-svar(ns)%Inh)*hh/24 - &
spar(ns)%PPb*svar(ns)%Pro*(24-hh)/24
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa*&
triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,atemp)* &
svar(ns)%Inh*hh/24
CASE(4)
! Quercus
! Promotor-Inhibitor model 12
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa* &
triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,atemp)* &
(1-svar(ns)%Inh)*hh/24
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa* &
triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,atemp)* &
svar(ns)%Inh*hh/24 + spar(ns)%PPb*(24-hh)/24
CASE(5, 11)
! Betula, Robinia
IF(spar(ns)%Phmodel==1) THEN
! Promotor-Inhibitor model 2
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa* &
triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,atemp)* &
(1-svar(ns)%Inh) - spar(ns)%PPb*svar(ns)%Pro*(24-hh)/24
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa* &
triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,atemp)*svar(ns)%Inh
END IF
END SELECT
Enddo
IF(spar(spec_new)%phmodel==4) THEN
svar(spec_new)%daybb = svar(spec_new)%ext_daybb
ELSE
svar(spec_new)%daybb = 181 + leapyear(time_cur)
ENDIF
endif
endif
endif
enddo
endif ! time
k = 0
do i=1,nspecies
if (helpanz(i) > 0) then
k = k + 1
anrspec = k
nrspec(k) = i
endif
ntr = svar(i)%sum_ntreea
if (svar(i)%sumNPP .gt. 1E-06) svar(i)%drIndAl = svar(i)%drIndAl / svar(i)%sumNPP
if (i .le. nspec_tree) then
IF(helpanz(i) > 0) THEN
anz_spec = anz_spec + 1
IF(helpdiam(i) < ntr) THEN
svar(i)%med_diam = SQRT(svar(i)%med_diam / (ntr - helpdiam(i)))
ENDIF
svar(i)%Ndem = svar(i)%Ndem / kpatchsize ! g per tree --> g/m2
svar(i)%Nupt = svar(i)%Nupt / kpatchsize ! g per tree --> g/m2
if (ntr .ne. 0) then
svar(i)%mean_height = svar(i)%mean_height / ntr
svar(i)%mean_diam = svar(i)%mean_diam / ntr
svar(i)%mean_jrb = svar(i)%mean_jrb / ntr
svar(i)%basal_area = pi*(ntr-helpdiam(i))*(svar(i)%med_diam*svar(i)%med_diam/40000)*10000/kpatchsize
else
svar(i)%sum_ntreea = 0.
endif
call calc_heidom_spec(i)
ENDIF
end if ! nspec_tree
! conversion kg/patch ---> kg/ha; N/patch ---> N/ha
helpntr = svar(i)%sum_nTreeA* 10000./kpatchsize
if(helpntr.eq.0 .and. svar(i)%sum_nTreeA.eq.1) then
svar(i)%sum_nTreeA = 1
else
svar(i)%sum_nTreeA = helpntr
end if
svar(i)%sum_bio = svar(i)%sum_bio * 10000./kpatchsize
svar(i)%fol = svar(i)%fol * 10000./kpatchsize
svar(i)%sap = svar(i)%sap* 10000./kpatchsize
svar(i)%hrt= svar(i)%hrt* 10000./kpatchsize
svar(i)%frt= svar(i)%frt* 10000./kpatchsize
svar(i)%totstem_m3= ( svar(i)%sap + svar(i)%hrt)/ (spar(i)%prhos*1000000) ! m3/ha
svar(i)%totsteminc = svar(i)%totsteminc * 10000./kpatchsize ! kg/ha
svar(i)%totsteminc_m3 = svar(i)%totsteminc_m3 * 10000./kpatchsize ! kg/ha
totsteminc_m3 = totsteminc_m3 + svar(i)%totsteminc_m3
totsteminc = totsteminc + svar(i)%totsteminc
end do
! new calculation of dominant height
call calc_heidom
if(anz_tree>0)then
if(lowtree<anz_tree) then
med_diam = sqrt(med_diam /(anz_tree-lowtree))
mean_diam = mean_diam /(anz_tree-lowtree)
mean_height = mean_height /(anz_tree-lowtree)
basal_area = pi*(anz_tree-lowtree)*(med_diam*med_diam/40000)*10000/kpatchsize
endif
else
if (hntr .ne. 0) then
med_diam = sqrt(med_diam /hntr)
mean_diam = mean_diam / hntr
mean_height = mean_height / hntr
else
med_diam = 0.
mean_diam = 0.
mean_height = 0.
endif
endif
end subroutine stand_bal_spec
!**************************************************************
subroutine class
use data_stand
use data_simul
use data_species
use data_par
implicit none
integer i,k
diam_class=0;height_class=0
diam_class_age=0.
diam_class_h = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
k = zeig%coh%species
if (k.ne.nspec_tree+2) then !exclusion of mistletoe
if(zeig%coh%diam<=dclass_w) then
diam_class(1,k)=diam_class(1,k)+zeig%coh%ntreea
diam_class_h(1,k) = diam_class_h(1,k) + zeig%coh%ntreea*zeig%coh%height
diam_class_age(1,k) = diam_class_age(1,k)+zeig%coh%x_age*zeig%coh%ntreea
end if
do i=2,num_class
if(zeig%coh%diam.le.(dclass_w + dclass_w*(i-1)) .and. zeig%coh%diam>(dclass_w + dclass_w*(i-2))) then
diam_class(i,k)=diam_class(i,k) + zeig%coh%ntreea
diam_class_h(i,k) = diam_class_h(i,k) + zeig%coh%ntreea*zeig%coh%height
diam_class_age(i,k) = diam_class_age(i,k)+zeig%coh%x_age*zeig%coh%ntreea
else if(zeig%coh%diam.gt. (dclass_w + dclass_w*(num_class-2))) then
diam_class(num_class,k)=diam_class(num_class,k) + zeig%coh%ntreea
diam_class_h(num_class,k) = diam_class_h(num_class,k) + zeig%coh%ntreea*zeig%coh%height
diam_class_age(num_class,k) = diam_class_age(num_class,k) + zeig%coh%x_age+zeig%coh%ntreea
exit
end if
enddo
if(zeig%coh%height<=100) height_class(1) = height_class(1)+zeig%coh%ntreea
if(zeig%coh%height>100.and.zeig%coh%height<500) height_class(2) = height_class(2)+zeig%coh%ntreea
do i=3,num_class-2
if(zeig%coh%height>(i+2)*100.and.zeig%coh%height<=(i+3)*100) height_class(i) = height_class(i)+zeig%coh%ntreea
enddo
if(zeig%coh%height>5000.and.zeig%coh%height<5500) height_class(num_class-1) = height_class(num_class-1)+zeig%coh%ntreea
if(zeig%coh%height>5500) height_class(num_class) = height_class(num_class)+zeig%coh%ntreea
endif!exclusion of mistletoe
zeig=>zeig%next
enddo
do i=1,num_class
do k=1,nspec_tree
if(diam_class(i,k).ne.0) diam_class_h(i,k) = (diam_class_h(i,k)/diam_class(i,k))*10000./kpatchsize
if(diam_class_age(i,k).ne.0.and.diam_class(i,k).ne.0 ) diam_class_age(i,k) =diam_class_age(i,k)/diam_class(i,k)
diam_class(i,k) = diam_class(i,k)*10000./kpatchsize
end do
end do
end subroutine class
!**************************************************************
subroutine classt
use data_stand
use data_simul
use data_species
implicit none
integer i,k
diam_class_t=0;height_class=0
diam_class_h = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
k = zeig%coh%species
if (k.ne.nspec_tree+2) then ! exclusion mistletoe
if(zeig%coh%diam<=dclass_w) then
diam_class_t(1,k)=diam_class_t(1,k)+zeig%coh%ntreed
end if
do i=2,num_class
if(zeig%coh%diam.le.(dclass_w + dclass_w*(i-1)) .and. zeig%coh%diam>(dclass_w + dclass_w*(i-2))) then
diam_class_t(i,k)=diam_class_t(i,k) + zeig%coh%ntreed
else if(zeig%coh%diam.gt. (dclass_w + dclass_w*(num_class-2))) then
diam_class_t(num_class,k)=diam_class_t(num_class,k) + zeig%coh%ntreed
exit
end if
enddo
endif !exclusion of mistletoe
zeig=>zeig%next
enddo
do i=1,num_class
do k=1,nspec_tree
diam_class_t(i,k)=diam_class_t(i,k)*10000./kpatchsize
end do
end do
end subroutine classt
!**************************************************************
subroutine class_man
use data_stand
use data_simul
use data_species
use data_manag
implicit none
integer i , k
real anz
diam_classm=0
diam_classm_h=0.
diam_class_mvol = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem.ne.0.or.(zeig%coh%ntreed.gt.0 .and. zeig%coh%diam.gt.tardiam_dstem)) then
if(zeig%coh%diam.le.tardiam_dstem) then
anz = zeig%coh%ntreem
else
anz = zeig%coh%ntreem + zeig%coh%ntreed
end if
k = zeig%coh%species
if(zeig%coh%diam<=dclass_w) then
diam_classm(1,k)=diam_classm(1,k)+anz
diam_classm_h(1,k) = diam_classm_h(1,k) + anz*zeig%coh%height
diam_class_mvol(1,k) = diam_class_mvol(1,k) +anz*(zeig%coh%x_sap + zeig%coh%x_hrt)
end if
if(zeig%coh%diam<=dclass_w*2.and.zeig%coh%diam.gt.dclass_w) then
diam_classm(2,k)=diam_classm(2,k)+anz
diam_classm_h(2,k) = diam_classm_h(2,k) + anz*zeig%coh%height
diam_class_mvol(2,k) = diam_class_mvol(2,k) + anz*(zeig%coh%x_sap + zeig%coh%x_hrt)
end if
do i=3,num_class
if(zeig%coh%diam.le.(dclass_w*2 + dclass_w*(i-2)) .and. zeig%coh%diam>(dclass_w*2 + dclass_w*(i-3))) then
diam_classm(i,k) = diam_classm(i,k) + anz
diam_classm_h(i,k) = diam_classm_h(i,k) + anz*zeig%coh%height
diam_class_mvol(i,k) = diam_class_mvol(i,k) + anz*(zeig%coh%x_sap + zeig%coh%x_hrt)
else if(zeig%coh%diam.gt. (dclass_w*2 + dclass_w*(num_class-3))) then
diam_classm(num_class,k)=diam_classm(num_class,k) + anz
diam_classm_h(num_class,k) = diam_classm_h(num_class,k) + anz*zeig%coh%height
diam_class_mvol(num_class,k) = diam_class_mvol(num_class,k) + anz*(zeig%coh%x_sap + zeig%coh%x_hrt)
end if
enddo
end if
zeig=>zeig%next
enddo
do i=1,num_class
do k=1,nspecies
if(diam_classm(i,k).ne.0) diam_classm_h(i,k) = diam_classm_h(i,k)/diam_classm(i,k)
if(diam_class_mvol(i,k).ne.0.) then
diam_class_mvol(i,k) = diam_class_mvol(i,k)/(spar(k)%prhos*1000000)*10000/kpatchsize
end if
diam_classm(i,k) = diam_classm(i,k)*10000./kpatchsize
end do
end do
end subroutine class_man
!**************************************************************
subroutine calc_heidom
use data_out
use data_simul
use data_stand
implicit none
real :: mh
integer :: nhelp, &
nh1,nh2, &
testflag=0, &
j
allocate (height_rank(anz_coh))
nh1=0
nh2=0
mh = 0
testflag = 0
nhelp = nint(kpatchsize/100)
if(anz_tree.le.nhelp) nhelp = anz_tree
! sorting by height of cohorts
call dimsort(anz_coh, 'height',height_rank)
if(anz_tree>1) then
do j=anz_coh, 1,-1
call dimsort(anz_coh, 'height',height_rank)
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ident.eq.height_rank(j)) then
nh2 = nh1
nh1 = nh1 + zeig%coh%ntreea
if(nh1.le. nhelp) then
mh = mh + zeig%coh%ntreea*zeig%coh%height
else
mh = mh + zeig%coh%height*( nhelp - nh2)
testflag=1
exit
end if
endif
zeig=>zeig%next
if(testflag.eq.1) exit
end do
if(testflag.eq.1) exit
if(nh1.eq.nhelp) exit
end do
if (nhelp.lt. nh1) then
hdom = mh/nhelp
else
hdom = mh/nh1
end if
end if
deallocate(height_rank)
end subroutine calc_heidom
!**************************************************************
subroutine calc_heidom_spec(ispec)
!*****************************************************
! species specific dominant height calculation
!*****************************************************
use data_out
use data_simul
use data_stand
implicit none
real :: mh
integer :: nhelp, &
nh1,nh2, &
testflag=0, &
j, &
ispec
allocate (height_rank(anz_coh))
hdom = 0
nh1=0
nh2=0
mh = 0
testflag = 0
! calculation of number of trees used for H100 ( 100/ ha = nhelp/ kpachtsize)
nhelp = nint(kpatchsize/100)
if(anz_tree.le.nhelp) nhelp = anz_tree
! sorting by height of cohorts
call dimsort(anz_coh, 'height',height_rank)
if(anz_tree>1) then
do j=anz_coh, 1,-1
call dimsort(anz_coh, 'height',height_rank)
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ident.eq.height_rank(j).and. zeig%coh%species.eq.ispec) then
nh2 = nh1
nh1 = nh1 + zeig%coh%ntreea
if(nh1.le. nhelp) then
mh = mh + zeig%coh%ntreea*zeig%coh%height
else
mh = mh + zeig%coh%height*( nhelp - nh2)
testflag=1
exit
end if
endif
zeig=>zeig%next
if(testflag.eq.1) exit
end do
if(testflag.eq.1) exit
if(nh1.eq.nhelp) exit
end do
if (nhelp.lt. nh1.and. nhelp.ne.0) then
hdom = mh/nhelp
else if(nh1.ne.0) then
hdom = mh/nh1
end if
else if(anz_tree.eq.1) then
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.ispec) hdom=zeig%coh%height
zeig=>zeig%next
end do
end if
deallocate(height_rank)
svar(ispec)%dom_height = hdom
end subroutine calc_heidom_spec
!**************************************************************
subroutine max_height(nrmax,anz,cohl)
use data_out
use data_simul
use data_stand
implicit none
integer :: nrmax
integer :: nrmax_h
integer :: anz, testflag,i
real :: help_h1, help_h2
integer,dimension(0:anz_coh) :: cohl
testflag=0
nrmax = -1
nrmax_h = -1
help_h2=0.
help_h1=0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
do i=0,anz-1
if(cohl(i).eq.zeig%coh%ident) then
testflag=1
endif
end do
if (testflag.eq.0) then
help_h2= zeig%coh%height
nrmax_h = zeig%coh%ident
if(help_h2.gt. help_h1) then
help_h1 = help_h2
nrmax = nrmax_h
end if
end if
zeig=>zeig%next
testflag = 0
end do
anz = anz +1
cohl(anz-1) = nrmax
end subroutine max_height
!**************************************************************
SUBROUTINE standup
! update of stand variables (LAI, cover, waldtyp)
USE data_par
USE data_stand
USE data_soil
USE data_species
use data_out
use data_simul
implicit none
integer i
REAL :: sumfol_can = 0.
REAL :: sumfol_sveg= 0.
REAL :: ntr, cover3
! estimating degree of covering
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time, ' standup'
cover3 = 0.
sumfol_can = 0.
sumfol_sveg= 0.
crown_area = 0.
do i = 1, anrspec
svar(nrspec(i))%crown_area = 0.
enddo
zeig=>pt%first
do
IF(.not.associated(zeig)) exit
if (zeig%coh%crown_area .ge. 0) then
ntr = zeig%coh%nTreeA
ns = zeig%coh%species
cover3 = cover3 + ntr * zeig%coh%crown_area
if (ns .le. nspec_tree) then
sumfol_can = sumfol_can + ntr * zeig%coh%x_fol
crown_area = crown_area + ntr * zeig%coh%crown_area
else
sumfol_sveg = sumfol_sveg + ntr * zeig%coh%x_fol
endif
svar(ns)%crown_area = svar(ns)%crown_area + ntr * zeig%coh%crown_area
endif
zeig=>zeig%next
end do
cover3 = cover3 / kpatchsize
anz_tree = 0
zeig=>pt%first
do
IF(.not.associated(zeig)) exit
ns=zeig%coh%species
if (ns .le. nspec_tree) then
zeig%coh%rel_fol = zeig%coh%ntreea * zeig%coh%x_fol/sumfol_can
ceppot_can = ceppot_can + zeig%coh%rel_fol * spar(ns)%ceppot_spec
anz_tree = anz_tree + zeig%coh%ntreea
else if (ns.eq.nspec_tree+1) then
zeig%coh%rel_fol = zeig%coh%ntreea * zeig%coh%x_fol/sumfol_sveg
ceppot_sveg = ceppot_sveg + zeig%coh%rel_fol * spar(ns)%ceppot_spec
endif
zeig=>zeig%next
end do
!Berechnung LAI und ceppot
ceppot_can = 0.
ceppot_sveg = 0.
LAI_can = 0.
LAI_sveg = 0.
DO i=1,anrspec
ns = nrspec(i)
IF (ns .le. nspec_tree) THEN
LAI_can = LAI_can + svar(ns)%act_sum_lai
ELSE
LAI_sveg = LAI_sveg + svar(ns)%act_sum_lai
ENDIF
ENDDO
DO i=1,anrspec
ns = nrspec(i)
IF (ns .le. nspec_tree) THEN
IF(LAI_can .gt. 0.) THEN
ceppot_can = ceppot_can + svar(ns)%act_sum_lai/LAI_can * spar(ns)%ceppot_spec
ELSE
ceppot_can = 0.
ENDIF
ELSE
IF(LAI_sveg .gt. 0.) THEN
ceppot_sveg = ceppot_sveg + svar(ns)%act_sum_lai/LAI_sveg * spar(ns)%ceppot_spec
ELSE
ceppot_sveg= 0.
ENDIF
END IF
ENDDO
if (LAI .gt. 1.) then
cover = cover3
else if (LAI .le. zero) then
cover = 0.1 * cover3
else
cover = LAI * cover3 ! to combine with leave surface; an Blattflaeche koppeln
endif
call wclas(waldtyp) ! forest type
END SUBROUTINE standup
!******************************************************************************
SUBROUTINE senescence
! update of senescence rates
USE data_stand
USE data_species
USE data_simul
IMPLICIT NONE
! senescence rates
zeig=>pt%first
DO
IF(.not.associated(zeig)) exit
if (zeig%coh%species.ne.nspec_tree+2) then ! exclude mistletoe from senescence
select case (flag_dis)
! case (1,2)
! zeig%coh%sfol = spar(zeig%coh%species)%psf * zeig%coh%x_fol + zeig%coh%x_fol_loss
! zeig%coh%sfrt = spar(zeig%coh%species)%psr * zeig%coh%x_frt + zeig%coh%x_frt_loss
case (0,1,2)
zeig%coh%sfol = spar(zeig%coh%species)%psf * zeig%coh%x_fol
zeig%coh%sfrt = spar(zeig%coh%species)%psr * zeig%coh%x_frt
end select
IF (zeig%coh%height.ge.thr_height .and.zeig%coh%species.LE. nspec_tree) THEN
zeig%coh%ssap = spar(zeig%coh%species)%pss * zeig%coh%x_sap
ELSE
zeig%coh%ssap = 0
if(zeig%coh%species.GT.nspec_tree) zeig%coh%ssap = spar(zeig%coh%species)%pss*zeig%coh%x_sap
ENDIF
end if !exclusion of mistletoe
zeig=>zeig%next
END DO
END SUBROUTINE senescence
!**************************************************************
SUBROUTINE litter
! Calculation of summation variables of litter fractions
use data_par
use data_out
use data_simul
use data_soil
use data_soil_cn
use data_species
use data_stand
implicit none
real hconvd
integer taxnr, i
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' litter'
zeig => pt%first
do while (associated(zeig))
taxnr = zeig%coh%species
if(taxnr.le.nspec_tree) then
totfol_lit_tree = totfol_lit_tree + zeig%coh%litC_fol
totfrt_lit_tree = totfrt_lit_tree + zeig%coh%litC_frt
end if
totfol_lit = totfol_lit + zeig%coh%litC_fol
totfrt_lit = totfrt_lit + zeig%coh%litC_frt
tottb_lit = tottb_lit + zeig%coh%litC_tb
totcrt_lit = totcrt_lit + zeig%coh%litC_crt
totstem_lit = totstem_lit + zeig%coh%litC_stem
zeig => zeig%next
enddo ! zeig (cohorts)
! litter biomass: x kg C/tree to kg/ha (n*x*1000g/(kPatchSize m2)/cpart==> kg/ha)
hconvd = (1000*gm2_in_kgha) / (kpatchsize * cpart) !
totfrt_lit = totfrt_lit * hconvd
totfol_lit = totfol_lit * hconvd
tottb_lit = tottb_lit * hconvd
totcrt_lit = totcrt_lit * hconvd
totstem_lit = totstem_lit * hconvd
totfol_lit_tree = totfol_lit_tree * hconvd
totfrt_lit_tree = totfrt_lit_tree * hconvd
do i = 1,nspec_tree
tottb_lit = tottb_lit + dead_wood(i)%C_tb(1)*gm2_in_kgha
totstem_lit = totstem_lit + dead_wood(i)%C_stem(1)*gm2_in_kgha
enddo
END subroutine litter
!**************************************************************
SUBROUTINE calc_ind_rep
USE data_stand
USE data_species
USE data_simul
implicit none
integer :: i
real :: hi
real, dimension(nspecies) :: rindex_spec
rindex1 = 0.
rindex2 = 0.
if(anz_spec.ne.0) then
hi = 1/real(anz_spec)
rindex_spec = 0.
do i = 1, nspecies
if(sumbio.ne.0) then
rindex_spec(i) = svar(i)%sum_bio/sumbio
end if
end do
rindex1 = 0.
rindex2 = 1.
do i = 1, nspecies
if(rindex_spec(i).ne.0) then
rindex1 = rindex1 + abs(hi -rindex_spec(i))
rindex2 = rindex2 * abs(hi -rindex_spec(i))
end if
end do
if(hi.ne.1) then
rindex1 = 1. - rindex1/(2*(1.-hi))
else
rindex1 = 0.
end if
rindex2 = rindex2**anz_spec
end if
END subroutine calc_ind_rep
!**************************************************************
SUBROUTINE overstorey
use data_out
USE data_stand
USE data_species
USE data_simul
implicit none
real,dimension(nspec_tree) :: mindbh, maxdbh, dminage, dmaxage
integer :: i, nrmin, taxnr, agedmin, agedmax
real :: dbhmin, dbhmax
integer :: anzoverst, nrmax
anzoverst = 0
mindbh=0.
do i =1,nspec_tree
call min_dbh(nrmin,dbhmin,agedmin, i)
mindbh(i) = dbhmin
dminage(i) = agedmin
call max_dbh(nrmax, dbhmax, agedmax, i)
maxdbh(i) = dbhmax
dmaxage(i) = agedmax
end do
if (time.eq.0) then
zeig=>pt%first
do
IF(.not.associated(zeig)) exit
taxnr = zeig%coh%species
if(taxnr .le.nspec_tree) then
if(zeig%coh%x_age.lt. (dminage(taxnr) +20) .and. dminage(taxnr).lt. dmaxage(taxnr)) then
zeig%coh%underst =2
end if
end if
zeig=>zeig%next
end do
else
zeig=>pt%first
do
IF(.not.associated(zeig)) exit
taxnr = zeig%coh%species
if(zeig%coh%height.gt. 130..and. zeig%coh%underst.eq.4) then
zeig%coh%underst = 2
end if
zeig=>zeig%next
end do
end if ! time
END SUBROUTINE overstorey
source_code/version_2.3_windows/stand_mort.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* - stand_mort *!
!* - int_mort intrinsic mortality rate *!
!* - stress_mort stress mortality rate *!
!* - int_mort_weib *!
!* *!
!* - Calculation of dead trees per cohort and species *!
!* deterministic approach *!
!* - relative mortality rate is determined by intr. mortality *!
!* and stress mortality *!
!* - stress mortality is calculated depending on *!
!* npp, ystress, yhealth *!
!* - intrinsic probability is optionally calculated on *!
!* age of cohort *!
!* - for each tree of the cohort mortality is decided *!
!* by means of the Mortality probability and *!
!* a uniformly distributed variable *!
!* *!
!* 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 *!
!* *!
!*****************************************************************!
! input variables:
! pro cohort NPP
! state variables:
! pro cohort nTreeA nTreeD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE stand_mort
USE data_stand
USE data_species
USE data_simul
USE data_manag
use data_out
use data_par
IMPLICIT NONE
!local variables
INTEGER :: flag_hgrowth
INTEGER :: taxnr
REAL :: intmort
REAL :: strmort
REAL :: totmort
REAL :: totmort_m
REAL :: besmort
REAL :: ntdead
REAL :: ntdead_m
REAL :: nhelp
REAL :: survpfunct
INTEGER :: hage
REAL :: besp1,besp2 ! parameters for besetting mortality
real :: help1, help2
real :: intmorth
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' stand_mort'
if (flag_standup .eq. 0) flag_standup = 1 ! call stand_balance later
ntdead=0.
nhelp=0.
ntdead_m=0
flag_hgrowth=0
sumvsdead = 0.
sumvsdead_m3 = 0.
besmort = 0
strmort = 0.
totmort = 0.
svar%sumvsdead = 0
svar%sumvsdead_m3 = 0.
zeig=>pt%first
DO
IF(.not.associated(zeig)) exit
IF (zeig%coh%height.ge.thr_height.and. zeig%coh%species.le.nspec_tree) then
taxnr=zeig%coh%species
IF(time.eq.1) then
zeig%coh%nta=zeig%coh%nTreeA
ELSE
IF (flag_mg.eq.1) then
IF (thin_year(act_thin_year-1).eq.(time-1)) zeig%coh%nta=zeig%coh%nTreeA
ELSE IF (flag_mg.eq.2) then
if(flag_adapm .eq. 1) zeig%coh%nta=zeig%coh%nTreeA
ENDIF
ENDIF
IF(zeig%coh%notViable) then
print*,time, zeig%coh%notViable
zeig%coh%nTreeD = zeig%coh%ntreeA
zeig%coh%nta = 0.
zeig%coh%ntreeA = 0
goto 1000
ENDIF
! calculation of stress and health indicator based on foliage biomass increment
hage = zeig%coh%x_age
IF (flag_hgrowth==0) THEN
IF(zeig%coh%fol_inc.le.0.0) then
zeig%coh%x_stress = zeig%coh%x_stress + 1
zeig%coh%x_health= 0
ELSE
zeig%coh%x_health = zeig%coh%x_health + 1
IF(zeig%coh%x_stress.eq.1.and. zeig%coh%x_health.gt.0) zeig%coh%nta = zeig%coh%ntreeA
IF(zeig%coh%x_stress .gt.0) zeig%coh%x_stress = zeig%coh%x_stress - 1
ENDIF
ENDIF
IF (flag_hgrowth==1) THEN
IF(zeig%coh%bio_inc.le.0.0) then
zeig%coh%x_stress = zeig%coh%x_stress + 1
zeig%coh%x_health= 0
ELSE
zeig%coh%x_health = zeig%coh%x_health + 1
IF(zeig%coh%x_stress.eq.1.and. zeig%coh%x_health.gt.0) zeig%coh%nta = zeig%coh%ntreeA
IF(zeig%coh%x_stress .gt.0) zeig%coh%x_stress = MAX(zeig%coh%x_stress - 3,0)
ENDIF
ENDIF
! calculation of relative mortality rates
! intrinsic mortality
! constant
call int_mort(taxnr,intmorth)
! age-depending using Weibull function
call int_mort_weib(taxnr,intmort,hage)
! stress mortality
IF(zeig%coh%x_stress.gt.0) then
IF(flag_hgrowth==0) strmort = weibal*spar(taxnr)%weibla*zeig%coh%x_stress**(weibal-1)
IF(flag_hgrowth==1) strmort = weibal*spar(taxnr)%weibla*(zeig%coh%x_stress/3.)**(weibal-1)
ELSE
strmort = 0.
ENDIF
!mortality depending on gross growth rate foliage
IF(strmort==0.0 .AND. flag_hgrowth==2) THEN
IF(zeig%coh%sfol/zeig%coh%gfol.GT.0.9) THEN
strmort=((zeig%coh%sfol/zeig%coh%gfol-0.9)*10.)**2
ELSE
ENDIF
ENDIF
if(strmort==0. .and. flag_hgrowth==3) then
help1 = 10**((log10(4.5)-log10((zeig%coh%x_sap + zeig%coh%x_hrt))*zeig%coh%ntreea)/1.5)
help2 = help1/zeig%coh%ntreea
end if
! mortality caused by besetting for oak
besp1= 0.018
besp2= 0.0216
if(zeig%coh%species.eq.4) then
if( zeig%coh%bes.le. 1.2) then
besmort = 0.
else
besmort = besp1*zeig%coh%bes- besp2
end if
else if (zeig%coh%species.eq.1.) then
if(zeig%coh%bes.le. 1.2) then
besmort = 0.
else
besmort = 0.04*zeig%coh%bes- 0.04
end if
end if
!total mortality rate depending on flag_mort
IF(flag_mort.eq.1) THEN
totmort = strmort
ELSE IF(zeig%coh%x_age.le.30) then
totmort=intmort+(1-intmort)*strmort
if(taxnr.eq.8) totmort = strmort
ELSE
totmort=intmort+(1-intmort)*strmort
ENDIF
! if species type oak then combination of stress mortality and besetting mortality
if(zeig%coh%species.eq.4.or.zeig%coh%species.eq.1) then
totmort = besmort + (1-besmort)* strmort
end if
survpfunct = exp(- spar(taxnr)%weibla * zeig%coh%x_stress**weibal)
ntdead = totmort*zeig%coh%nTreeA
IF(totmort.GT.1) CALL error_mess(time,"totmort greater 1: ",totmort)
! calculation of real stem number
zeig%coh%nta = zeig%coh%nta - ntdead
! calculation of integer stem number
zeig%coh%nTreeD = zeig%coh%nTreeA-NINT(zeig%coh%nta)
zeig%coh%nTreeA = NINT(zeig%coh%nta)
IF(zeig%coh%nTreeA.lt.1.) zeig%coh%nTreeA=0.
if (zeig%coh%mistletoe.eq.1) then ! in case Mist.infect. tree dies, number of mistletoes dies, too
totmort_m = zeig%coh%nTreeD/(zeig%coh%nTreeD+zeig%coh%nTreeA) ! share of trees removed of total trees. used for the share of mistletoe that dies
ntdead_m = 1. !flag
end if
! calculation of the litter pool of a died tree of a cohort
1000 IF (zeig%coh%ntreeD.ne.0) then
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreeD*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreeD*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreeD*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreeD*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreeD*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreeD*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreeD*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreeD*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
if(flag_mg.ne.0) then
! if(thin_dead.eq.0.and.thin_flag1(1).lt.0.) then
if(thin_dead.eq.0) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
end if
else if(zeig%coh%diam.le.tardiam_dstem.or. flag_mg.eq.0) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
else if(zeig%coh%diam.gt.tardiam_dstem.and.flag_mg.ne.0.or.flag_mg.eq.5) then
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
else if(thin_dead.eq.1.and.zeig%coh%diam.le.tardiam_dstem) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
end if
ENDIF
ENDIF
zeig=>zeig%next
ENDDO
! if tree cohort with mistletoe changed, change number of mistletoes too
if (ntdead_m.eq.1.) then
zeig => pt%first
do while (associated(zeig))
if (zeig%coh%species.eq.nspec_tree+2) then
zeig%coh%nta=zeig%coh%nTreeA
ntdead_m = totmort_m*zeig%coh%nTreeA
zeig%coh%nta = zeig%coh%nta - ntdead_m
zeig%coh%nTreeD = zeig%coh%nTreeA-NINT(zeig%coh%nta)
zeig%coh%nTreeA = NINT(zeig%coh%nta)
IF(zeig%coh%nTreeA.lt.1.) then
zeig%coh%nTreeA=0.
flag_mistle=0 !set flag mistletoe back to zero
ENDIF
endif
zeig=>zeig%next
enddo ! zeig (cohorts)
end if
ntdead_m=0.
! recalculation sumvsdead
sumvsdead = sumvsdead * 10000./kpatchsize ! kg/patch ---> ! kg/ha
sumvsdead_m3 = sumvsdead_m3 * 10000./kpatchsize ! kg/patch ---> ! kg/ha
cumsumvsdead = cumsumvsdead + sumvsdead
END SUBROUTINE stand_mort
SUBROUTINE int_mort(taxnr,intmort)
USE data_species
IMPLICIT NONE
REAL :: intmort
INTEGER :: taxnr
intmort=1.-exp(-spar(taxnr)%intr)
END SUBROUTINE int_mort
SUBROUTINE int_mort_weib(taxnr,intmort,hage)
USE data_species
USE data_stand
USE data_simul
IMPLICIT NONE
REAL :: intmort, weibla_int
INTEGER :: taxnr
INTEGER :: hage
! Weibull functions depending on age
weibla_int = -log(0.01)/(spar(taxnr)%max_age**weibal_int)
intmort = weibal_int*weibla_int*(hage)**(weibal_int-1.)
END SUBROUTINE int_mort_weib
source_code/version_2.3_windows/stand_regen.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* subroutine for regeneration *!
!* including the SR: *!
!* - gener_seed *!
!* - seed_ini *!
!* - simseed *!
!* - growth_seed *!
!* - mort_seed *!
!* *!
!* 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 stand_regen
USE data_simul
USE data_stand
IMPLICIT NONE
flag_standup = 2
CALL mort_seed
CALL gener_seed
END SUBROUTINE stand_regen
SUBROUTINE gener_seed
USE data_stand
USE data_species
USE data_simul
use data_out
USE data_plant
USE data_soil
IMPLICIT NONE
real :: seedla ! leaf area of all seedling cohorts
real :: laiseed ! lai ----"------
integer :: nseed ! number of generated seeds
real :: redseed
integer :: i
integer, dimension(5) :: agemin, seedpot
real, dimension(5,3) :: latg
real :: help,help1, help2
real :: pequal
integer :: hlayer
integer :: flag_reg_help
TYPE(coh_obj), POINTER :: p
DATA latg /1.,0.3,0.1,1.,0.1,1.,0.9,0.5,1.,0.5,1.,1.,0.9,1.,0.9/
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' gener_seed'
flag_reg_help = 0
agemin = 0
seedpot = 0
seedla = 0.
help1 = 0.
! calculation of leafarea of all seedling cohorts
SELECT CASE (flag_reg)
! according to FORGRA (Vincent Kint)
CASE (30)
call random_number(pequal)
DO i= 1, nspecies
nseed = 0
p =>pt%first
DO WHILE (ASSOCIATED(p))
if(p%coh%species.eq.i) then
if (i.eq.1) then
agemin(i) = 50 + (30* pequal)
call random_number(pequal)
seedpot(i) = 810*pequal
else if(i.eq.3) then
agemin(i) = 15 + (45* pequal)
call random_number(pequal)
seedpot(i) = 1000*pequal
else if(i.eq.4) then
agemin(i) = 15 +(35* pequal)
call random_number(pequal)
seedpot(i) = 1125*pequal
else if(i.eq.5) then
agemin(i) = 10 +(10* pequal)
call random_number(pequal)
seedpot(i) = 8750*pequal
end if
if(p%coh%x_age.ge. agemin(i).and.p%coh%diam.gt.0.) then
nseed = nseed + seedpot(i)*(p%coh%ntreem + p%coh%ntreea)
end if
end if
p => p%next
END DO ! cohort
help2 = irelpool_ll
if(help2.lt.0) help2 =0
if(help2.eq. 0.) then
redseed = 0.
else if( help2.gt.0. .and. help2.le.latg(i,1)) then
redseed = help2*latg(i,1)/0.4
else if ( help2.gt.latg(i,1).and. help2.le.latg(i,2)) then
redseed = help2*latg(i,2)/0.6
else if ( help2.gt.latg(i,2).and. help2.le.latg(i,3)) then
redseed = help2*latg(i,3)/0.8
else if(help2.gt.latg(i,3)) then
redseed = help2* latg(i,3)
end if
nseed = nseed * redseed
! for birch 1 year old saplings are used
if (i.eq.5) then
numplant(i) = nseed
flag_reg = 14
if(nseed.ne.0) call planting
flag_reg= 0
else
call seed_multi(nseed,i)
end if
END DO ! species
CASE(1,2,3)
p =>pt%first
DO WHILE (ASSOCIATED(p))
if(p%coh%height.lt.thr_height) then
seedla = seedla + p%coh%t_leaf*p%coh%ntreea
help1 = help1 + p%coh%x_fol*p%coh%ntreea
end if
p => p%next
END DO
! calculation LAI of lowest_layer
laiseed=seedla/kpatchsize
DO i=1,nspecies
IF (spar(i)%regflag.eq.1) THEN
CALL simseed(i,nseed)
IF(laiseed.lt.1) THEN
! reduction of seedling number nseed depending on light module and free space in the lowest_layer
SELECT CASE (flag_light)
CASE(1)
IF(flag_reg.ne.3) THEN
CALL seed_ini(nseed,i)
ELSE
CALL seed_multi(nseed,i)
END IF
CASE (2)
if (anz_coh.eq. 0) then
if(time.eq.1) then
hlayer = 0
else
hlayer = 1
end if
else
hlayer = lowest_layer -1
end if
help = vstruct(hlayer)%Irel
if (help.lt.0.05) help = 0
IF(help.eq.0) THEN
nseed = 0
ELSE
nseed = nseed*help
IF(flag_reg.ne.3) THEN
CALL seed_ini(nseed,i)
ELSE
CALL seed_multi(nseed,i)
END IF
END IF
CASE(3)
redseed= bgpool_ll
nseed = nseed*redseed
IF(flag_reg.ne.3) THEN
CALL seed_ini(nseed,i)
ELSE
CALL seed_multi(nseed,i)
END IF
CASE(4)
if(i.gt.5) then
redseed= irelpool_ll
else
! according to FORGRA, not for all species (i=1,5)
help2 = irelpool_ll
if(help2.lt. 0.01) then
redseed = 0.
else if( help2.gt.0.01 .and. help2.le.latg(i,1)) then
redseed = help2*latg(i,1)/0.4
else if ( help2.gt.latg(i,1).and. help2.le.latg(i,2)) then
redseed = help2*latg(i,2)/0.6
else if ( help2.gt.latg(i,2).and. help2.le.latg(i,3)) then
redseed = help2*latg(i,3)/0.8
else if(help2.gt.latg(i,3)) then
redseed = help2* latg(i,3)
end if
end if
nseed = redseed * nseed
IF(flag_reg.ne.3) THEN
CALL seed_ini(nseed,i)
ELSE
if (i.eq.5) then
numplant(i) = nseed
flag_reg_help = flag_reg
flag_reg = 14
if(nseed.ne.0) call planting
flag_reg = flag_reg_help
else
CALL seed_multi(nseed,i)
end if
END IF
END SELECT
ELSE
nseed = 0.
END IF
ELSE
nseed=0.
END IF
END DO
END SELECT ! flag_reg
END subroutine gener_seed
SUBROUTINE simseed(specnum,nseed)
USE data_species
use data_simul
use data_stand
IMPLICIT NONE
REAL :: pequal
INTEGER :: nseed,specnum
REAL :: seedmax
! calculation of max. seedrate of patch from max. seedrate per m2
seedmax=spar(specnum)%seedrate*kpatchsize
CALL random_number(pequal)
CALL random_number(pequal)
nseed=-seedmax*alog(1.-pequal)
IF (flag_mg ==4 .and. time.eq.1) THEN
nseed = NINT(spar(specnum)%seedrate*kpatchsize)
ELSE IF(flag_mg ==4.and. time.gt.1)THEN
nseed = 0
END IF
end
SUBROUTINE seed_ini(nseed,nsp)
USE data_species
use data_stand
use data_help
use data_out
use data_simul
use data_soil
IMPLICIT NONE
integer :: nseed, nr, j
integer :: nsp
REAL :: shoot
REAL :: x1,x2,xacc,root
REAL :: rtflsp
REAL :: troot2
TYPE(cohort) ::tree_ini
external weight
external rtflsp
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' seed_ini'
IF(nseed.eq.0) RETURN
hnspec = nsp
max_coh = max_coh + 1
! nullify of all elements
call coh_initial (tree_ini)
tree_ini%ident = max_coh
tree_ini%species = nsp
tree_ini%ntreea = nseed
tree_ini%nta = nseed
tree_ini%x_age = 1
mschelp = spar(nsp)%seedmass/1000. ! g ---> kg
x1 = 0.
x2 = 0.1
xacc = (1.0e-10) * (x1+x2)/2
root = rtflsp(weight,x1,x2,xacc)
tree_ini%x_sap = root
shoot = root*1000. ! [kg]
tree_ini%x_fol= (spar(nsp)%seeda*(tree_ini%x_sap** spar(nsp)%seedb)) ![kg]
tree_ini%x_frt = tree_ini%x_fol ! [kg]
tree_ini%med_sla = spar(nsp)%psla_min + spar(nsp)%psla_a*0.5
tree_ini%t_leaf = tree_ini%med_sla* tree_ini%x_fol ! [m-2]
tree_ini%ca_ini = tree_ini%t_leaf
! initialize pheno state variables
IF(spar(tree_ini%species)%Phmodel==1) THEN
tree_ini%P=0
tree_ini%I=1
ELSE
tree_ini%P=0
tree_ini%I=0
tree_ini%Tcrit=0
END IF
! tranformation of shoot biomass kg --> mg
if(nsp.ne.2)tree_ini%height = spar(nsp)%pheight1*(shoot*1000.)**spar(nsp)%pheight2 ! [cm] calculated from shoot biomass (mg); berechnet aus shoot biomass (mg)
if(nsp.eq.2) tree_ini%height = 10**(spar(nsp)%pheight1+ spar(nsp)%pheight2*LOG10(shoot*1000.)+ &
spar(nsp)%pheight3*(LOG10(shoot*1000.))**2)
IF(nseed.ne.0.) then
IF (.not. associated(pt%first)) THEN
ALLOCATE (pt%first)
pt%first%coh = tree_ini
NULLIFY(pt%first%next)
! root distribution
call root_depth (1, pt%first%coh%species, pt%first%coh%x_age, pt%first%coh%height, pt%first%coh%x_frt, pt%first%coh%x_crt, nr, troot2, pt%first%coh%x_rdpt, pt%first%coh%nroot)
pt%first%coh%nroot = nr
do j=1,nr
pt%first%coh%rooteff = 1. ! assumption for the first use
enddo
do j=nr+1, nlay
pt%first%coh%rooteff = 0. ! layers with no roots
enddo
ELSE
ALLOCATE(zeig)
zeig%coh = tree_ini
zeig%next => pt%first
pt%first => zeig
! root distribution
call root_depth (1, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot)
zeig%coh%nroot = nr
do j=1,nr
zeig%coh%rooteff = 1. ! assumption for the first use
enddo
do j=nr+1, nlay
zeig%coh%rooteff = 0. ! layers with no roots
enddo
END IF
anz_coh=anz_coh+1
END IF
END SUBROUTINE seed_ini
SUBROUTINE growth_seed
USE data_stand
USE data_species
USE data_simul
USE data_par
use data_out
IMPLICIT NONE
REAL :: lambdaf = 0., & ! partitioning functions
lambdas = 0., &
lambdar = 0., &
NPP = 0., & ! annual NPP
F = 0., & ! state variables: foliage,
S = 0., & ! shoot biomass,
R = 0., & ! fine roots,
H = 0., & ! total tree height
FNew, SNew, & ! new states
RNew, &
sigmaf = 0., & ! current leaf activity rate
sigman = 0., & ! current root activity rate
betar = 0., &
ar = 0
REAL :: Sf, & ! senescence rates
Sr, &
Gf, & ! growth rates
Gs, &
Gr
real :: pab,helpdr,helpsum
TYPE(coh_obj), POINTER :: p
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' growth_seed'
flag_standup = 2 ! call stand_balance and root_distribution later
p=>pt%first
DO
if(.not.associated(p)) exit
if( p%coh%height.lt.thr_height.and. p%coh%species.le.nspec_tree) then
ns = p%coh%species
F = p%coh%x_fol
S = p%coh%x_sap
R = p%coh%x_frt
NPP = p%coh%NPP
IF (flag_reg .eq. 2) NPP = p%coh%NPPpool ! [kg]
H = p%coh%height
Sf = p%coh%sfol
Sr = p%coh%sfrt
! only allocate if enough NPP is available
1 IF (NPP>1.0E-9.or. NPP.ge.(Sf+Sr).and.(sr+Sf)>1.0E-9) THEN
! calculate leaf activity based on net PS and leaf mass
sigmaf = NPP/F
! calculate root activity based on drought index
helpdr= p%coh%drIndAl / p%coh%nDaysGr
IF (flag_sign.eq.1) THEN
sigman = amax1(spar(ns)%sigman*10*(((5.-spar(ns)%stol)*1.-p%coh%crown_area)/(5.-spar(ns)%stol)*1.),spar(ns)%sigman) * p%coh%drIndAl / p%coh%nDaysGr
ELSE
sigman = spar(ns)%sigman * p%coh%drIndAl / p%coh%nDaysGr
END IF
! auxiliary variables for fine roots
if(helpdr.lt.0.001) ar = 0.
ar = spar(ns)%pcnr * sigmaf / sigman
betar = (Sr - R + ar*(F-Sf)) / NPP
! calculate coefficients for roots and foliage and shoot
select case (ns)
case (1)
pab = 0.487
case(2)
pab = 0.826
case(3)
pab=1.9
case(4)
pab=1.002
! Pinus contorta
case(6)
! Gholz
pab = 0.236
! Populus tremula
case(8)
pab = 0.3233
case(9)
! Pinus halepensis
pab = 1.42335
case(10)
! pseudotsuga menziesii
pab = 0.8515
case(11)
! Robinia
pab = 0.8594
end select
lambdaf = (pab*(1-betar)+ (Sf/NPP))/(1 + pab*(1. + ar))
lambdar = ar * lambdaf + betar
lambdas = 1.- lambdaf - lambdar
! consistency
ELSE
lambdaf = 0.
lambdas = 0.
lambdar = 0.
END IF
if ( lambdas.lt.0.) then
lambdas = 0.
lambdaf = (1.-betar)/(ar+1)
lambdar = 1.-lambdaf
if( lambdar.lt.0) then
lambdar=0
lambdaf=1
end if
if(lambdaf.lt.0) then
lambdaf =1
lambdar = 0.
end if
endif
helpsum = lambdaf + lambdar+ lambdas
Gf = lambdaf*NPP
Gr = lambdar*NPP
Gs = lambdas*NPP
p%coh%gfol = Gf
p%coh%gfrt = Gr
p%coh%gsap = Gs
! update of state vector
FNew = F + Gf - Sf
SNew = S + Gs
RNew = R + Gr - Sr
p%coh%x_fol = FNew
p%coh%x_sap = SNew
p%coh%x_frt = RNew
p%coh%fol_inc_old = p%coh%fol_inc
p%coh%fol_inc = Gf - Sf
p%coh%stem_inc = Gs
! update height and shoot base diameter (regression functions from Schall 1998)
IF(ns.ne.2) p%coh%height = spar(ns)%pheight1* (snew*1000000.) **spar(ns)%pheight2
IF(ns.eq.2) p%coh%height = 10**(spar(ns)%pheight1+ spar(ns)%pheight2*LOG10(snew*1000000.)+ &
spar(ns)%pheight3*(LOG10(snew*1000000.))**2)
p%coh%height_ini = p%coh%height
! update foliage area, parameter med_sla
SELECT CASE (flag_light)
CASE (1:2)
p%coh%med_sla = spar(ns)%psla_min + spar(ns)%psla_a*(1.- vstruct(lowest_layer)%irel)
CASE(3,4)
p%coh%med_sla = spar(ns)%psla_min + spar(ns)%psla_a*(1.-irelpool(lowest_layer))
END SELECT
! total leaf area of a tree in this cohort [m**2]
p%coh%ca_ini = p%coh%med_sla * p%coh%x_fol
! update age -now not necessary this is done in stand_bal
p%coh%notViable = (FNew <= 0.) .OR. (SNew <= 0.) .OR. &
(RNew <= 0.)
p%coh%litC_fol = p%coh%litC_fol + p%coh%ntreea * Sf * cpart
p%coh%litC_frt = p%coh%litC_frt + p%coh%ntreea * Sr * cpart
! with species specific N content and reallocation factor (see species.par)
p%coh%litN_fol = p%coh%litN_fol + p%coh%ntreea * Sf * cpart * spar(ns)%reallo_fol / spar(ns)%cnr_fol
p%coh%litN_frt = p%coh%litN_frt + p%coh%ntreea * Sr * cpart * spar(ns)%reallo_frt / spar(ns)%cnr_frt
end if ! seedling cohort test
p=> p%next
END DO
END SUBROUTINE growth_seed
SUBROUTINE mort_seed
USE data_species
USE data_simul
use data_stand
use data_par
use data_out
IMPLICIT NONE
integer :: taxnr
integer :: hage
real :: intmort
real :: strmort
real :: totmort
real :: ntdead
real :: ntahelp
TYPE(coh_obj), POINTER :: p
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' mort_seed'
p=>pt%first
DO
IF(.not.associated(p)) EXIT
IF(p%coh%height.lt.thr_height) THEN
IF(p%coh%notViable) then
PRINT*,time, p%coh%notViable
p%coh%ntreed = p%coh%ntreea
p%coh%ntreea = 0
ENDIF
END IF
p => p%next
END DO
p => pt%first
DO
IF(.not.associated(p)) EXIT
IF(p%coh%height.lt.thr_height .and. p%coh%species.le.nspec_tree) THEN
taxnr = p%coh%species
if(p%coh%ntreea .eq.0) goto 1000
hage = p%coh%x_age
IF( p%coh%fol_inc.le.0.) THEN
p%coh%x_stress = p%coh%x_stress +1
p%coh%x_health = 0
ELSE
p%coh%x_health = p%coh%x_health +1
IF(p%coh%x_stress .gt. 0.) p%coh%x_stress = p%coh%x_stress -1
ENDIF
! intrinsic mortality
CALL int_mort_weib(taxnr, intmort, hage)
! stress mortality
IF(p%coh%x_stress.gt.0) THEN
strmort = weibal*spar(taxnr)%weibla*p%coh%x_stress**(weibal-1)
ELSE
strmort = 0.
ENDIF
totmort=intmort+(1-intmort)*strmort
! calculation of real number of dying stems
ntdead = totmort*p%coh%ntreeA
! update of real stem number nta and number of dead stems
p%coh%nta = p%coh%nta -ntdead
p%coh%nTreeD = p%coh%nTreeA-NINT(p%coh%nta)
! help variable for comparison
ntahelp = p%coh%nTreeA
! update of integer stem number
p%coh%nTreeA = NINT(p%coh%nta)
! update of integer stem number
if(p%coh%nta.lt.1.) p%coh%nTreeA=0.
! update of real stem number if integer stem number was changed
if (ntahelp .ne. p%coh%nTreeA ) p%coh%nta = p%coh%nTreeA
1000 if (p%coh%ntreeD.ne.0) then
p%coh%litC_fol = p%coh%litC_fol + p%coh%ntreeD*(1.-spar(taxnr)%psf)*p%coh%x_fol*cpart
p%coh%litN_fol = p%coh%litN_fol + p%coh%ntreeD*((1.-spar(taxnr)%psf)*p%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
p%coh%litC_frt = p%coh%litC_frt + p%coh%ntreeD*p%coh%x_frt*cpart
p%coh%litN_frt = p%coh%litN_frt + p%coh%ntreeD*p%coh%x_frt*cpart/spar(taxnr)%cnr_frt
p%coh%litC_tb = p%coh%litC_tb + p%coh%ntreeD*p%coh%x_tb*cpart
p%coh%litN_tb = p%coh%litN_tb + p%coh%ntreeD*p%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
p%coh%litC_crt = p%coh%litC_crt + p%coh%ntreeD*p%coh%x_crt*cpart
p%coh%litN_crt = p%coh%litN_crt + p%coh%ntreeD*p%coh%x_crt*cpart/spar(taxnr)%cnr_tbc
p%coh%litC_stem = p%coh%litC_stem + p%coh%ntreeD*(p%coh%x_sap)*cpart
p%coh%litN_stem = p%coh%litC_stem/spar(taxnr)%cnr_stem
endif
END IF
p => p%next
ENDDO
END SUBROUTINE mort_seed
source_code/version_2.3_windows/statistik.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* *!
!* statistical analysis of model quality *!
!* *!
!* Author: F. Suckow *!
!* *!
!* contains: *!
!* residuen *!
!* statistik *!
!* mean (n, arr) *!
!* variance (n, meanv, arr) *!
!* correl (n, meanv1, arr1, meanv2, arr2) *!
!* sumsq (n, arr) *!
!* stat_mon *!
!* *!
!* 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 residuen (ip)
use data_mess
implicit none
integer i,j, ires, ip
! calculate and save residues, with date, simulation as well as measurement value
! Residuen berechnen, mit Datum, Sim.- und Messwert speichern
if (ip .eq. 1) then
allocate (val(imkind))
do i=1,imkind
ires = 0
val(i)%tkind = tkind
allocate (val(i)%day(1:anz_val))
allocate (val(i)%year(1:anz_val))
allocate (val(i)%resid(1:anz_val))
allocate (val(i)%sim(1:anz_val))
allocate (val(i)%mess(1:anz_val))
val(i)%day = -99
val(i)%year = -99
val(i)%resid = -9999.0
val(i)%mkind = sim_kind(i)
do j = 1,anz_val
if (mess2(j,i) .gt. -9000.0 .and. sim1(j,i) .gt. -9000.0) then
ires = ires + 1
val(i)%day(ires) = stz(1,j)
val(i)%year(ires) = stz(2,j)
val(i)%resid(ires)= sim1(j,i) - mess2(j,i)
val(i)%sim(ires) = sim1(j,i)
val(i)%mess(ires) = mess2(j,i)
else
endif
enddo
val(i)%imes = ires
enddo
else
do i=1,imkind
ires = 0
val(i)%resid = -9999.9
val(i)%sim = -9999.9
do j = 1,anz_val
if (mess2(j,i) .gt. -9000.0 .and. sim1(j,i) .gt. -9000.0) then
ires = ires + 1
val(i)%resid(ires)= sim1(j,i) - mess2(j,i)
val(i)%sim(ires) = sim1(j,i)
else
endif
enddo
enddo
endif
END SUBROUTINE residuen
!**************************************************************
SUBROUTINE statistik
use data_mess
use data_simul
implicit none
integer imt ! aktueller Messwert-Typ
real, external :: mean, variance, correl, sumsq
integer i, n, nhelp
real help, h1, h2
real, allocatable, dimension(:):: arr, arrs, arrm, harr
real:: avs, & ! mean value simulation; Mittelwert Simulation
mins, & ! Minimum Simulation
maxs, & ! Maximum Simulation
stdevs, & ! standard deviation simulation; Standardabweichung Simulation
varis, & ! scattering of simulation; Streuung Simulation
varcos, & ! coefficient of variation for simulation; Variationskoeffizient Simulation
avm, & ! mean value measurements; Mittelwert Messwerte
minm, & ! minimum value measurements; Minimum Messwerte
maxm, & ! maximum value measurements; Maximum Messwerte
stdevm, & ! standard deviation measurements; Standardabweichung Messwerte
varim, & ! scattering of measurements; Streuung Messwerte
varcom, & ! coefficient of variation of measurements; Variationskoeffizient Messwerte
corrco, & ! coefficient of correlation; Korrelationskoeffizient
rsq, & ! coefficient of determination; Bestimmtheitsmass
avr, & ! mean error residues; Mittlerer Fehler Residuen
minr, & ! minimum residues; Minimum Residuen
maxr, & ! maximum residues; Maximum Residuen
stdevr, & ! standard deviation residues; Standardabweichung Residuen
varir, & ! scattering of residues; Streuung Residuen
varcor, & ! coefficient of variation residues; Variationskoeffizient Residuen
nme, & ! normalised mean error; Normalisierter mittlerer Fehler
mae, & ! mean absolute error of residues; Mittlerer absoluter Fehler Residuen
nmae, & ! normalised mean absolute error; Normalisierter mittlerer absoluter Fehler
sse , & ! sum of squared errors; Fehlerquadratsumme
rmse, & ! Root mean square error
nrmse, & ! Normalised root mean square error
pme, & ! mean procental error; Mittlerer prozentualer Fehler
prmse, & ! mean squared procentual error; Mittlerer quadratischer prozentualer Fehler
tic, & ! Theilsch imbalance coefficient; Theilscher Ungleichheitskoeffizient
meff ! Model efficiency (Medlyn et al. 2005)
do imt = 1, imkind
n = val(imt)%imes
if (n .gt. 0) then
allocate (arr(n))
allocate (arrs(n))
allocate (arrm(n))
allocate (harr(n))
! simulation
arrs = val(imt)%sim(1:n)
avs = mean(n, arrs)
mins = minval(arrs)
maxs = maxval(arrs)
varis = variance(n, avs, arrs)
if (varis .ge. 0.) then
stdevs = sqrt(varis)
else
stdevs = 0.
endif
if (avs .ne. 0.) then
varcos = stdevs / avs
else
varcos = -9999.0
endif
! observed
arrm = val(imt)%mess(1:n)
avm = mean(n, arrm)
minm = minval(arrm)
maxm = maxval(arrm)
varim = variance(n, avm, arrm)
if (varim .ge. 0.) then
stdevm = sqrt(varim)
else
stdevm = 0.
endif
! residuals
arr = val(imt)%resid(1:n)
avr = mean(n, arr)
minr = minval(arr)
maxr = maxval(arr)
varir = variance(n, avr, arr)
if (varir .ge. 0.) then
stdevr = sqrt(varir)
else
stdevr = 0.
endif
if (avr .ne. 0.) then
varcor = stdevr / avr
else
varcor = -9999.0
endif
corrco = correl(n, avs, arrs, avm, arrm)
if (corrco .ge. -1.) then
rsq = corrco * corrco
rsq_av = rsq_av + rsq
else
imk_rsq = imk_rsq - 1
rsq = -9999.0
endif
if (avs .ne. 0.) then
nme = (avm - avs) / avs
nme_av = nme_av + nme
else
imk_nme = imk_nme - 1
nme = -9999.0
endif
mae = mean(n, abs(arr))
sse = sumsq(n, arr)
rmse = sqrt(sse / n)
if (avm .ne. 0.) then
varcom = stdevm / avm
nrmse = rmse / abs(avm)
nrmse_av = nrmse_av + nrmse
nmae = mae / abs(avm)
nmae_av = nmae_av + nmae
else
imk_nrmse = imk_nrmse - 1
imk_nmae = imk_nmae - 1
varcom = -9999.0
nrmse = -9999.0
nmae = -9999.0
endif
nhelp = n
do i = 1, n
if (arrm(i) .ne. 0.) then
harr(i) = abs(arr(i)/arrm(i))
else
nhelp = nhelp -1
harr(i) = 0
endif
enddo
pme = mean(nhelp, harr)
prmse = sumsq(nhelp, harr)
prmse = sqrt(prmse / nhelp)
tic = sse / (sumsq(n, arrs) + sumsq(n, arrm))
h1=sumsq(n, arr)
harr = arrm-avm
h2=sumsq(n, harr)
if (h2.ne.0) then
meff = 1. - (h1 / h2)
else
meff=1
end if
! Mittelwert
pme_av = pme_av + pme
prmse_av = prmse_av + prmse
tic_av = tic_av + tic
meff_av = meff_av + meff
deallocate (arr)
deallocate (arrm)
deallocate (arrs)
deallocate (harr)
write (unit_stat, '(I5,2X, A20,1X,A10,I8,1X,30E13.5)') ip, site_name(ip), val(imt)%mkind, val(imt)%imes, &
avr, minr, maxr, stdevr, varir, varcor, nme, mae, nmae, sse, rmse, nrmse, pme, prmse, tic,meff, corrco, rsq, &
avs, mins, maxs, stdevs, varis, varcos, avm, minm, maxm, stdevm, varim, varcom
endif
enddo
END SUBROUTINE statistik
!**************************************************************
REAL FUNCTION mean (n, arr)
integer n, i
real, dimension(n):: arr
real help
help = 0.
do i = 1,n
help = help + arr(i)
enddo
mean = help / n
END FUNCTION mean
!**************************************************************
REAL FUNCTION variance (n, meanv, arr)
integer n, i
real, dimension(n):: arr
real meanv, help, xx
help = 0.
if (n .gt. 1) then
do i = 1,n
xx = arr(i) - meanv
help = help + xx * xx
enddo
variance = help / (n -1)
else
variance = -9999.0
endif
END FUNCTION variance
!**************************************************************
REAL FUNCTION correl (n, meanv1, arr1, meanv2, arr2)
integer n, i
real, dimension(n):: arr1, arr2
real meanv1, meanv2, help1, help2, help3, xx1, xx2
help1 = 0.
help2 = 0.
help3 = 0.
do i = 1,n
xx1 = arr1(i) - meanv1
xx2 = arr2(i) - meanv2
help1 = help1 + xx1 * xx2
help2 = help2 + xx1 * xx1
help3 = help3 + xx2 * xx2
enddo
if ((help2 .gt. 1.E-06) .and. (help3 .gt. 1.E-06)) then
correl = help1 / sqrt(help2*help3)
else
correl = -9999.0
endif
END FUNCTION correl
!**************************************************************
REAL FUNCTION sumsq (n, arr)
integer n, i
real, dimension(n):: arr
real help
help = 0.
do i = 1,n
help = help + arr(i) * arr(i)
enddo
sumsq = help
END FUNCTION sumsq
!**************************************************************
Subroutine stat_mon
! Statistics of monthly values, derived from daily observed values
use data_mess
use data_out
use data_simul
implicit none
integer i, j, k, l
integer dd, mm, yy, doy, yanz, arranz
integer :: outunit ! output unit
character(250) text, filename
character(20) idtext, datei, vunit, obskind
character(150) htext
real, allocatable, dimension(:):: helparr ! help array with montly values of one month for all years
real, allocatable, dimension(:,:):: help_mon ! array with monthly values for each year, year
real, allocatable, dimension(:,:):: help_day ! array with mean daily values per month for each year, year
integer, allocatable, dimension(:,:):: help_num ! array with number of measurement values for each year, year
yanz = mtz(2,imess) - mtz(2,1) + 1
if (.not. allocated(unit_mon)) then
allocate(unit_mon(imkind))
allocate(unit_mon_stat(imkind))
allocate(helparr(yanz))
endif
if (.not. allocated(help_mon)) then
allocate(help_mon(12,yanz))
allocate(help_day(12,yanz))
allocate(help_num(12,yanz))
endif
do k = 1, imkind
help_mon = 0.0
help_num = 0
obskind = sim_kind(k)
filename = trim(dirout)//trim(site_name(ip))//'_'//trim(obskind)//'_mon_obs'//'.out'
unit_mon(k) = getunit()
open(unit_mon(k),file=filename,status='replace')
! Calculate mmonthly sums
do j = 1, imess
doy = mtz(1,j)
yy = mtz(2,j)
call TZINDA(dd,mm,yy,doy)
yy = mtz(2,j) - mtz(2,1) + 1
if (mess1(j,k) .Gt. -9990.) then
if (sim_kind(k) .eq. 'AET') then
if (mess1(j,k) .lt. 0.) then
mess_info = '# negative AET set to zero'
mess1(j,k) = 0. ! avoid negative AET
endif
endif
help_mon(mm,yy) = help_mon(mm,yy) + mess1(j,k)
help_num(mm,yy) = help_num(mm,yy) + 1
endif
enddo ! j
do j = 1, yanz
do i = 1,12
if (help_num(i,j) .gt. 0) then
help_day(i,j) = help_mon(i,j) / help_num(i,j)
else
help_mon(i,j) = -9999.
help_day(i,j) = -9999.
endif
enddo
enddo
! Output of monthly sums
select case (trim(obskind))
case ('AET')
vunit = 'mm'
case ('GPP', 'NPP', 'TER')
vunit = 'g C/m'
case ('Snow')
return
case default
vunit = ' '
end select
write (unit_mon(k), '(A)') '# Monthly sum, daily mean of month and number of values per month of observed '//trim(obskind)
write (unit_mon(k), '(A)') '# '//trim(vunit)
write (unit_mon(k), '(A)', advance='no') '# Year'
do i = 1,12
write (unit_mon(k), '(A8,I2)', advance='no') trim(obskind)//'_',i
enddo
write (unit_mon(k), '(A)')
l = 0
do j = mtz(2,1), mtz(2,imess)
l = l + 1
write (unit_mon(k), '(A,I6,12F10.2)') 'sum ', j, (help_mon(i,l), i=1,12)
write (unit_mon(k), '(A,I6,12F10.2)') 'daily mean ', j, (help_day(i,l), i=1,12)
write (unit_mon(k), '(A,I6,12I10)') 'number ', j, (help_num(i,l), i=1,12)
enddo
! Statistics
filename = trim(dirout)//trim(site_name(ip))//'_'//trim(obskind)//'_mon_obs_stat'//'.res'
outunit = getunit()
open(outunit,file=filename,status='replace')
write (outunit, '(A)') '# Statistics over all years for each monthly sum and daily mean per month of '//trim(obskind)
write (outunit, '(A)') '# '//trim(vunit)
write (outunit, '(A, I6)') '# Simulation period (years): ', year
write (outunit, '(A)') '# site_id Month number Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
write (outunit, '(A)') 'monthly sum'
do i = 1,12
arranz = 0
do j = 1,yanz
if (help_mon(i,j) .gt. -9990.) then
arranz = arranz + 1
helparr(arranz) = help_mon(i,j)
endif
enddo ! j
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! nur letzte 20 Zeichen schreiben
write (outunit, '(A20,I5,I8)', advance = 'no') idtext, i, arranz
call calc_stat(arranz, helparr, outunit)
enddo ! i
write (outunit, '(A)') ' '
write (outunit, '(A)') 'daily mean per month'
do i = 1,12
arranz = 0
do j = 1,yanz
if (help_day(i,j) .gt. -9990.) then
arranz = arranz + 1
helparr(arranz) = help_day(i,j)
endif
enddo ! j
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! nur letzte 20 Zeichen schreiben
write (outunit, '(A20,I5,I8)', advance = 'no') idtext, i, arranz
call calc_stat(arranz, helparr, outunit)
enddo ! i
enddo ! k
continue
End Subroutine stat_mon
\ No newline at end of file
source_code/version_2.3_windows/target_thin.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutine *!
!* target thinning - *!
!* thinning routine with given values of biomass per *!
!* thinning year as target values *!
!* targetm i given in kg DW/ha *!
!* *!
!* 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 target_thinning(i)
use data_stand
use data_manag
use data_simul
use data_species
use data_par
implicit none
real :: targetm ! target value of stem biomass
real :: dbhmin=0, &
dbhmin_us = 0, &
wpa=0, & ! Weibull parameter
wpa_us , &
wpb=0, & ! -"-
wpb_us, &
wpc=0, & ! -"-
d63=0, &
d63_us, &
help=0, &
pequal, &
tdbh=0, &
bas_area=0, &
bas_help=0., &
rtarget_help=0, &
target_help1=0,&
dbh_h =0, &
db_l = 0., &
db_u = 0., &
d_est=0., &
w_kb=0., &
stembiom, &
stembiom_us = 0. , &
stembiom_re = 0. , &
stembiom_all = 0. , &
diff, &
mdiam, &
mdiam_us
integer :: nrmin, &
lowtree, &
undertree, &
flagth, &
taxnr, &
counth, &
min_id, &
max_id, &
ih1,ih2,ncoh, &
coun1
! auxilary for thinning routine 4: selective thinning
integer :: count,i, &
idum , third, ipot, isel, ih
integer,dimension(0:anz_coh) :: cohl
integer, dimension(anz_coh) :: id_pot
real :: h1, h2
real,external :: gasdev
real:: ran0
! reacalculation of target to kg DW/patch
h1 = 0.
h2 = 0.
count = 0
cohl = -1
flagth = 0
coun1 = 0
help=0.
lowtree=0
undertree = 0
anz_tree_dbh = 0
bas_area = 0.
! stem biomass of overstorey
stembiom = 0.
! stem biomass of understorey
stembiom_us = 0.
stembiom_all = 0.
if (time.eq.73.and. ip .eq.87) then
stembiom = 0
end if
taxnr = thin_spec(i)
mdiam = 0.
mdiam_us = 0.
! calculation of mean diameter (correspondung to med_diam) and basal area of stand
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Modification for V Kint: no test for diameter
IF((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.0) THEN
! overstorey
stembiom = stembiom + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
help = help + zeig%coh%ntreea*(zeig%coh%diam**2)
bas_area = bas_area + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
if( zeig%coh%diam>0) then
anz_tree_dbh = anz_tree_dbh + zeig%coh%ntreea
mdiam = mdiam + zeig%coh%ntreea * (zeig%coh%diam**2)
end if
! Trees with DBH=0 for population and per species; Baeume mit DBH =0 fuer Bestand und pro Spezie
ELSE IF( (zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) THEN
! seedings/regeneration
stembiom_re = stembiom_re + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
lowtree = lowtree + zeig%coh%ntreea
ELSE if((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) THEN
! understorey
stembiom_us = stembiom_us + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
mdiam_us = mdiam_us + zeig%coh%ntreea * (zeig%coh%diam**2)
undertree = undertree + zeig%coh%ntreea
ENDIF
zeig => zeig%next
ENDDO
! mean diamteer for over and understorey
stembiom_all = stembiom + stembiom_us
if(anz_tree_dbh.ne.0) mdiam = sqrt(mdiam/real(anz_tree_dbh))
if(undertree.ne.0) mdiam_us = sqrt(mdiam_us/undertree)
third = nint(anz_tree_dbh*0.333333)
anz_tree_ha = nint(anz_tree_dbh*10000./kpatchsize)
IF(anz_tree>0)THEN
if(lowtree<anz_tree) help = sqrt(help/(anz_tree-lowtree))
ENDIF
! setting of aux. variable target_help
rtarget_help = stembiom_all
! tending
if(thin_tysp(i).eq.4.or.(stembiom_re.ne.0. .and. stembiom_all.eq.0)) then
rtarget_help = stembiom_re
end if
! Umrechnung in Biomasse pro patch in kg
targetm = target_mass(i)*1000*kpatchsize/10000.
! target value of biomass
if(thin_tysp(i).eq.4 .or.(stembiom_re.ne.0. .and. stembiom_all.eq.0) ) then
! tending
targetm = stembiom_re - targetm*stembiom_re
else
end if
if( targetm.eq.1) targetm = 0.
! targetm (kg DW/patch)
! cuttting
if (targetm.eq.0.)then
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
!tending of regeneration
else if(thin_tysp(i).eq.4) then
min_id = 1000
max_id = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1) then
ih1 = zeig%coh%ident
if(ih1.lt.min_id) min_id = ih1
ih2 = zeig%coh%ident
if (ih2.gt.max_id) max_id = ih2
end if
zeig=> zeig%next
end do
target_help1 = 0.
do
call random_number(pequal)
ncoh = min_id +(max_id-min_id)*pequal
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1.and. zeig%coh%ident.eq.ncoh ) then
zeig%coh%ntreea = zeig%coh%ntreea - 1
zeig%coh%nta = zeig%coh%nta-1
zeig%coh%ntreem = zeig%coh%ntreem +1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig=>zeig%next
end do
diff = targetm - rtarget_help
if(diff.lt.0.01) exit
end do
else IF ( targetm .ne. 0.) then
if(target_mass(i).lt.1.) then
targetm = target_mass(i) * rtarget_help
end if
! different thinnings from below and above
select case(thin_tysp(i))
case(1)
! moderate lower thinning;
d_est = 1.02
w_kb = 1.8
case(2)
! intensive lower thinning;
d_est = 1.03
w_kb = 1.5
case(3)
! high thinning;
d_est = 1.04
w_kb = 1.2
end select
! calculation of Weibull-Parameter
call min_dbh_overs(nrmin,dbhmin,taxnr)
call min_dbh_unders(nrmin,dbhmin_us, taxnr)
bas_help = bas_area
wpa = dbhmin
wpa_us = dbhmin_us
d63 = mdiam*d_est
d63_us = mdiam_us * d_est
wpb = (d63 - wpa)/ w_kb
wpb_us = (d63_us-wpa_us)/w_kb
wpc = 2
if (thin_tysp(i).eq. 3) then
! starting with overstorey!, continuing with the understorey
if(targetm.lt.(stembiom_all-stembiom)) then
! total removing of overstorey
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
! defining the remaining part of stem biomass which has to be remove
rtarget_help = stembiom_us
! understorey
!selection of trees for thinning
if(mdiam_us.ne.0) then
! start understorey thinning
do
call random_number(pequal)
tdbh = wpa_us + wpb_us*(-log(1.-pequal))**(1./wpc)
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt.10000) exit
end if
zeig=> zeig%next
END DO ! list of potential thinned tees cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
count = count +1
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(count.gt.100000) exit
end do ! understorey thinning
end if ! mediam_us.ne.0
else ! targetm.lt.(stembiom_all-stembiom)
!selection of trees for thinning
do
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.2*dbh_h
db_u = dbh_h + 0.2*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
call random_number(pequal)
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
! write(1234,*) time, ipot, pequal, isel
! if(isel.eq.0) isel =1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(coun1.gt.100000) exit
end do ! total thinning
end if !targetm.lt.(stembiom_all-stembiom)
end if ! thintype 3
if(thin_tysp(i).eq.1.or.thin_tysp(i).eq.2) then
if(targetm.lt.(stembiom_all-stembiom_us)) then
! total removing of understorey
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.2) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
! definging the remaining part of stem biomass which has to remove
rtarget_help = stembiom
if(mdiam.ne.0) then
! additional thinning from the overstorey
counth = 0
do
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.3*dbh_h
db_u = dbh_h + 0.3*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(counth.gt.100000) exit
end do ! overstorey thinning
end if ! mdiam.ne.0
else ! targtem.lt.(stembiom_all-stembiom_us)
! first thinning from the understorey
!selection of trees for thinning
if(mdiam_us.ne.0) then
do
call random_number(pequal)
tdbh = wpa_us + wpb_us*(-log(1.-pequal))**(1./wpc)
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth .gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1 .or. (stembiom_all-stembiom_us).eq.rtarget_help) exit
if(coun1.gt.100000) exit
end do ! understorey thinning
end if ! mdiam_us
end if
end if !! thin_tysp.eq.1 or. thin-tysp.eq.2
END IF ! all thinnings and tending
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem>0.and.zeig%coh%species.eq.taxnr) then
! all parts without stems of trees are input for litter
h1 = zeig%coh%ntreem*zeig%coh%x_fol
h2 = h2 + h1
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreem*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreem*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
endif
zeig=>zeig%next
enddo
END SUBROUTINE target_thinning
source_code/version_2.3_windows/target_thin_bas.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutine *!
!* target thinning - *!
!* thinning routine with given values of stem number per *!
!* thinning year as target values *!
!* rtargetm i given inN/ha *!
! or target value 0 < rtargetm <1 for *!
!* basal area reduction *!
!* *!
!* 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 target_thinning_bas(i)
use data_stand
use data_manag
use data_simul
use data_species
use data_par
implicit none
real :: rtargetm=0. ! target value of stem biomass
integer :: target_help = 0.
real :: dbhmin=0, &
dbhmin_us = 0, &
wpa=0, & ! Weibull parameter
wpb=0, & ! -"-
wpc=0, & ! -"-
d63=0, &
help=0, &
pequal, &
tdbh=0, &
bas_area=0, &
bas_help=0., &
bas_area_us = 0., &
bas_area_test=0., &
target_help1=0, &
dbh_h =0, &
db_l = 0., &
db_u = 0., &
d_est=0., &
w_kb=0., &
stembiom, &
stembiom_us, &
stembiom_re, &
stembiom_all, &
diff, &
mdiam, &
mdiam_us
integer :: nrmin, &
nrmin_us, &
lowtree, &
undertree, &
flagth, &
taxnr, &
counth, &
min_id, &
max_id, &
ih1,ih2,ncoh, &
coun1
integer :: flag_bas
! auxilarity for thinning routine 4: selective thinning
integer :: count, i, &
idum,third, ipot, isel, ih
integer,dimension(0:anz_coh) :: cohl
integer, dimension(anz_coh) :: id_pot
real :: h1, h2 , tar_h
real,external :: gasdev
real:: ran0
! reacalculation of target to kg DW/patch
flag_bas=0
h1 = 0.
h2 = 0.
count = 0
cohl = -1
flagth = 0
coun1 = 0
help=0.
lowtree=0
undertree = 0
anz_tree_dbh = 0
bas_area = 0.
bas_area_us = 0.
! stem biomass of overstorey
stembiom = 0.
! stem biomass of understorey
stembiom_us = 0.
stembiom_all = 0.
tar_h = 300.
if (time.eq.73.and. ip .eq.87) then
stembiom = 0
end if
taxnr = thin_spec(i)
mdiam = 0.
mdiam_us = 0.
! calculation of mean diameter (correspondung to med_diam) and basal area of stand
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Modification for V Kint: no test for diameter
IF((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.0) THEN
! overstorey
stembiom = stembiom + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
help = help + zeig%coh%ntreea*(zeig%coh%diam**2)
bas_area = bas_area + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
if( zeig%coh%diam>0) then
anz_tree_dbh = anz_tree_dbh + zeig%coh%ntreea
mdiam = mdiam + zeig%coh%ntreea * (zeig%coh%diam**2)
end if
! Trees with DBH=0 for populations and per species; Baeume mit DBH =0 fuer Bestand und pro Spezie
ELSE IF( (zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) THEN
! seedings/regeneration
stembiom_re = stembiom_re + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
lowtree = lowtree + zeig%coh%ntreea
ELSE if((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2.and.zeig%coh%underst.eq.2) THEN
! understorey
stembiom_us = stembiom_us + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
mdiam_us = mdiam_us + zeig%coh%ntreea * (zeig%coh%diam**2)
undertree = undertree + zeig%coh%ntreea
bas_area_us = bas_area_us + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
ENDIF
zeig => zeig%next
ENDDO
! mean diameter for over and understorey
stembiom_all = stembiom + stembiom_us
if(anz_tree_dbh.ne.0) mdiam = sqrt(mdiam/real(anz_tree_dbh))
if(undertree.ne.0) mdiam_us = sqrt(mdiam_us/undertree)
third = nint(anz_tree_dbh*0.333333)
anz_tree_ha = nint(anz_tree_dbh*10000./kpatchsize)
anz_tree = anz_tree_dbh + undertree
IF(anz_tree>0)THEN
if(lowtree<anz_tree) help = sqrt(help/(anz_tree-lowtree))
ENDIF
! new basal area management
if(target_mass(i).le.1) then
flag_bas=1
if(thin_stor(i).eq.0) then
target_help = bas_area
else
target_help = bas_area_us
end if
else
if(thin_stor(i).eq.0) then
target_help = anz_tree_dbh
else
target_help = undertree
end if
end if
! tending
if(thin_tysp(i).eq.4) then
target_help = bas_area_us
end if
if(target_mass(i).gt.1) then
rtargetm = target_mass(i)*kpatchsize/10000.
else
rtargetm = target_mass(i)
end if
! cuttting
if (rtargetm.eq.0.)then
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
!tending of regeneration
else if(thin_tysp(i).eq.4) then
rtargetm = target_mass(i) * target_help
min_id = 1000
max_id = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1 .or. zeig%coh%underst.eq.2) then
ih1 = zeig%coh%ident
if(ih1.lt.min_id) min_id = ih1
ih2 = zeig%coh%ident
if (ih2.gt.max_id) max_id = ih2
end if
zeig=> zeig%next
end do
target_help1 = 0.
do
call random_number(pequal)
ncoh = min_id +(max_id-min_id)*pequal
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1 .or.zeig%coh%underst.eq.2 .and. zeig%coh%ident.eq.ncoh ) then
zeig%coh%ntreea = zeig%coh%ntreea - 1
zeig%coh%nta = zeig%coh%nta-1
zeig%coh%ntreem = zeig%coh%ntreem +1
if (flag_bas.eq.1) then
target_help=target_help - zeig%coh%diam*zeig%coh%diam*pi/4
else
target_help = target_help - zeig%coh%ntreea
end if
exit
end if
zeig=>zeig%next
end do
diff = target_help - rtargetm
if(diff.lt.0.01) exit
end do
! different thinnings from below and above
else IF ( rtargetm .ne. 0.) then
if(target_mass(i).lt.1.) then
rtargetm = target_mass(i) * target_help
! No management if rtargetm=1
else if (rtargetm.eq.1) then
return
endif
select case(thin_tysp(i))
case(1)
! medium lower thinning
d_est = 1.02
w_kb = 2.5
case(2)
! strong lower thinning
d_est = 1.03
w_kb = 1.5
case(3)
! High thinning
d_est = 1.04
w_kb = 1.2
end select
! calculation of Weibull-Parameter
call min_dbh_overs(nrmin,dbhmin,taxnr)
call min_dbh_unders(nrmin_us,dbhmin_us, taxnr)
! write (7777,*) time, nrmin,nrmin_us, dbhmin, dbhmin_us
bas_help = bas_area
if (thin_stor(i).eq.0) then
wpa = dbhmin
else
wpa = dbhmin_us
end if
if (thin_stor(i).eq.0) then
d63 = mdiam*d_est
else
d63 = mdiam_us * d_est
end if
wpb = (d63 - wpa)/ w_kb
wpc = 2
wpc = 0.8
if ((thin_tysp(i).ne.4) .and. rtargetm.lt.target_help) then
!selection of trees for thinning
DO ! 11 begin selecting
IF (thin_stor(i) .eq. 2 .and. nrmin_us .eq. -1) EXIT ! exit thinning loop if there are no trees with dbh in understory an understory shall be thinned
coun1=coun1+1
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO ! 2 list preparing
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%notviable.eqv. .TRUE.) then
if(flag_mort.eq.0) then
id_pot(ipot)=zeig%coh%ident
ipot=ipot + 1
endif
else if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! 2 list of potential thinned tree cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
call random_number(pequal)
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO ! 3 removing one tree
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
if(zeig%coh%notviable.eqv..TRUE.) then
if(flag_mort.eq.0) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea=0
zeig%coh%nta=0
if (flag_bas.eq.1) then
target_help = target_help - (zeig%coh%diam**2)*pi/4
else
target_help = target_help -1
endif
end if
else
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
if(flag_bas.eq.1) then
! write (8888,*) 'target_help', time, target_help, rtargetm
target_help = target_help - (zeig%coh%diam**2)*pi/4
else
target_help = target_help -1
end if
exit
! end if
end if
end if
zeig =>zeig%next
ENDDO ! 3 thinning of one tree
end if
diff = target_help - rtargetm
if(diff.le.0.1) exit
if(coun1.gt.100000) exit
ENDDO ! 11 total thinning
end if ! thintype 1,2,3
END IF ! all thinnings and tending
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreea>0.and.zeig%coh%species.eq.taxnr) then
bas_area_test = bas_area_test + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
end if
zeig=>zeig%next
end do
! adding biomasses to litter pools depending on stage of stand
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem>0.and.zeig%coh%species.eq.taxnr) then
! all parts without stems of trees are input for litter
h1 = zeig%coh%ntreem*zeig%coh%x_fol
h2 = h2 + h1
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreem*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreem*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
endif
zeig=>zeig%next
enddo
call class_man
END SUBROUTINE target_thinning_bas
source_code/version_2.3_windows/target_thinstemnum.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutine *!
!* target thinning - *!
!* thinning routine with given values of biomass per *!
!* thinning year as target values *!
!* rtargetm i given in kg DW/ha *!
!* *!
!* 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 target_thinning_OC(i)
use data_stand
use data_manag
use data_simul
use data_species
use data_par
implicit none
real :: rtargetm=0. ! target value of stem biomass
integer :: target_help = 0.
real :: dbhmin=0, &
dbhmin_us = 0, &
wpa=0, & ! Weibull parameter
wpb=0, & ! -"-
wpc=0, & ! -"-
d63=0, &
help=0, &
pequal, &
tdbh=0, &
bas_area=0, &
bas_help=0., &
target_help1=0,&
dbh_h =0, &
db_l = 0., &
db_u = 0., &
d_est=0., &
w_kb=0., &
stembiom, &
stembiom_us, &
stembiom_re, &
stembiom_all, &
diff, &
mdiam, &
mdiam_us
integer :: nrmin, &
nrmin_us, &
lowtree, &
undertree, &
flagth, &
taxnr, &
counth, &
min_id, &
max_id, &
ih1,ih2,ncoh, &
coun1
! auxilarity for thinning routine 4: selective thinning
integer :: count, i, &
idum,third, ipot, isel, ih
integer,dimension(0:anz_coh) :: cohl
integer, dimension(anz_coh) :: id_pot
real :: h1, h2 , tar_h
real,external :: gasdev
real:: ran0
! reacalculation of target to kg DW/patch
h1 = 0.
h2 = 0.
count = 0
cohl = -1
flagth = 0
coun1 = 0
help=0.
lowtree=0
undertree = 0
anz_tree_dbh = 0
bas_area = 0.
! stem biomass of overstorey
stembiom = 0.
! stem biomass of understorey
stembiom_us = 0.
stembiom_all = 0.
tar_h = 300.
if (time.eq.73.and. ip .eq.87) then
stembiom = 0
end if
taxnr = thin_spec(i)
mdiam = 0.
mdiam_us = 0.
! calculation of mean diameter (correspondung to med_diam) and basal area of stand
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Modification for V Kint: no test for diameter
IF((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.0) THEN
! overstorey
stembiom = stembiom + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
help = help + zeig%coh%ntreea*(zeig%coh%diam**2)
bas_area = bas_area + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
if( zeig%coh%diam>0) then
anz_tree_dbh = anz_tree_dbh + zeig%coh%ntreea
mdiam = mdiam + zeig%coh%ntreea * (zeig%coh%diam**2)
end if
! Trees with DBH=0 for populations and per species; Baeume mit DBH =0 fuer Bestand und pro Spezie
ELSE IF( (zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) THEN
! seedings/regeneration
stembiom_re = stembiom_re + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
lowtree = lowtree + zeig%coh%ntreea
ELSE if((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) THEN
! understorey
stembiom_us = stembiom_us + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
mdiam_us = mdiam_us + zeig%coh%ntreea * (zeig%coh%diam**2)
undertree = undertree + zeig%coh%ntreea
ENDIF
zeig => zeig%next
ENDDO
! mean diameter for over and understorey
stembiom_all = stembiom + stembiom_us
if(anz_tree_dbh.ne.0) mdiam = sqrt(mdiam/real(anz_tree_dbh))
if(undertree.ne.0) mdiam_us = sqrt(mdiam_us/undertree)
third = nint(anz_tree_dbh*0.333333)
anz_tree_ha = nint(anz_tree_dbh*10000./kpatchsize)
anz_tree = anz_tree_dbh + undertree
IF(anz_tree>0)THEN
if(lowtree<anz_tree) help = sqrt(help/(anz_tree-lowtree))
ENDIF
if(thin_stor(i).eq.0) then
target_help = anz_tree_dbh
else
target_help = undertree
end if
! tending
if(thin_tysp(i).eq.4) target_help = stembiom_re
if(target_mass(i).gt.1) then
rtargetm = target_mass(i)*kpatchsize/10000.
else
rtargetm = target_mass(i)
end if
! target value of biomass
if(thin_tysp(i).eq.4) then
rtargetm = stembiom_re - rtargetm*stembiom_re
else
end if
! cuttting
if (rtargetm.eq.0.)then
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
!tending of regeneration
else if(thin_tysp(i).eq.4) then
min_id = 1000
max_id = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1) then
ih1 = zeig%coh%ident
if(ih1.lt.min_id) min_id = ih1
ih2 = zeig%coh%ident
if (ih2.gt.max_id) max_id = ih2
end if
zeig=> zeig%next
end do
target_help1 = 0.
do
call random_number(pequal)
ncoh = min_id +(max_id-min_id)*pequal
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1.and. zeig%coh%ident.eq.ncoh ) then
zeig%coh%ntreea = zeig%coh%ntreea - 1
zeig%coh%nta = zeig%coh%nta-1
zeig%coh%ntreem = zeig%coh%ntreem +1
target_help = target_help - zeig%coh%ntreea
exit
end if
zeig=>zeig%next
end do
diff = rtargetm -target_help
if(diff.lt.0.01) exit
end do
! different thinnings from below and above
else IF ( rtargetm .ne. 0.) then
if(target_mass(i).lt.1.) then
rtargetm = target_mass(i) * target_help
! No management if rtargetm=1
else if (rtargetm.eq.1) then
return
endif
select case(thin_tysp(i))
case(1)
! medium lower thinning
d_est = 1.02
w_kb = 2.5
case(2)
! strong lower thinning
d_est = 1.03
w_kb = 1.5
case(3)
! High thinning
d_est = 1.04
w_kb = 1.2
end select
! calculation of Weibull-Parameter
call min_dbh_overs(nrmin,dbhmin,taxnr)
call min_dbh_unders(nrmin_us,dbhmin_us, taxnr)
write (7777,*) time, nrmin, nrmin_us, dbhmin, dbhmin_us
bas_help = bas_area
if (thin_stor(i).eq.0) then
wpa = dbhmin
else
wpa = dbhmin_us
end if
if (thin_stor(i).eq.0) then
d63 = mdiam*d_est
else
d63 = mdiam_us * d_est
end if
wpb = (d63 - wpa)/ w_kb
wpc = 2
wpc = 0.8
if ((thin_tysp(i).ne.4) .and. rtargetm.lt.target_help) then
!selection of trees for thinning
do
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%notviable.eqv. .TRUE.) then
if(flag_mort.eq.0) then
id_pot(ipot)=zeig%coh%ident
ipot=ipot + 1
endif
else if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential thinned tree cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
call random_number(pequal)
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
if(zeig%coh%notviable.eqv..TRUE.) then
if(flag_mort.eq.0) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea=0
zeig%coh%nta=0
coun1=coun1+1
target_help = target_help - zeig%coh%ntreem
endif
else
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
target_help = target_help - 1
exit
end if
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = target_help - rtargetm
if(diff.le.0.1) exit
if(coun1.gt.100000) exit
end do ! total thinning
end if ! thintype 1,2,3
END IF ! all thinnings and tending
! adding biomasses to litter pools depending on stage of stand
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem>0.and.zeig%coh%species.eq.taxnr) then
! all parts without stems of trees are input for litter
h1 = zeig%coh%ntreem*zeig%coh%x_fol
h2 = h2 + h1
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreem*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreem*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
endif
zeig=>zeig%next
enddo
call class_man
END SUBROUTINE target_thinning_OC