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
Showing
with 68 additions and 632 deletions
File added
No preview for this file type
No preview for this file type
No preview for this file type
File deleted
No preview for this file type
No preview for this file type
File added
File added
File added
......@@ -121,7 +121,7 @@ module data_soil
! arrays of given root distribution (defined input)
real, allocatable, save, dimension(:) :: root_fr ! root fraction per soil layer
! dp_rfr ! depth of root fraction / cm
! yearly fine root loss after Rasse et al. 2001
integer :: rdepth_kind ! kind of calculation of root depth
real, allocatable, dimension(:) :: wat_left ! auxiliary variable for coh%watleft to determin annual sum of available water in soil layer boardering on root zone
......@@ -298,15 +298,15 @@ module data_soil_t
! Variables and parameters for soil temperature calculation
integer flag_surf ! calculation of soil surface temperature
! 0 - old version
! 1 - new ersion with explicit surface temperature
integer :: flag_surf = 0 ! calculation of soil surface temperature
! 0 - surface temperature equals temperature of first layer
! 1 - with explicit surface temperature
real temps_surf ! soil surface temperature
real hflux_surf ! soil heat flux at soil surface
real temps_surf ! soil surface temperature
real hflux_surf ! soil heat flux at soil surface
! model parameters
real :: C0 = 0.76, & ! Faltungskoeff.
real :: C0 = 0.76, & ! coefficients for calculation of surface temperature
C1 = 0.05, &
C2 = 0.3
......
......@@ -50,11 +50,12 @@ if (hum .le. 0.) then
else if (hum .gt. 100.) then
hum = 100.
endif
if (press .gt. 0.) then
if (prs(i,j) .gt. 0.) then
press = prs(i,j)
else
press = 1013.
endif
rad = rd(i,j)
wind = wd(i,j)
if (wind .lt. 0.) wind = 0.5
......
......@@ -118,7 +118,8 @@ SUBROUTINE INITIA
! end of declaration section
!******************************************************************************
ncl1 = 60
!ncl1 = 60
ncl1=40
allocate (zheigh(ncl1), zbhd(ncl1), zhbc(ncl1), nz(ncl1))
allocate (smaldc(ncl1), bigdc(ncl1))
print *,' '
......@@ -133,7 +134,7 @@ WRITE(*,'(A)',advance='no') ' ***Make your choice: '
READ *, data_flag
print *,' '
clwdth=2 !set diameter class-class width
clwdth=15 !set diameter class-class width
corr_la=1. !standard value for leaf area correction in stands of high sum of crown projection areas
mixed_tot_ca=0. !sum of crown projection area for mixed stands
pass = 1 !counter for number of passes through calculation loop for mixed stands
......@@ -256,7 +257,7 @@ CASE(1)
IF (datasets=='multi') THEN
select_lines=.false.
fl_num=0
if(infile=='input/hyyti_ini_0616.txt') then
ALLOCATE(ngroups(10000))
numstand= 0
......@@ -291,7 +292,6 @@ if(infile=='input/hyyti_ini_0616.txt') then
iF(baum(i).EQ.22) ngroups(nlines)%taxid=6 ! Larix
iF(baum(i).EQ.23) ngroups(nlines)%taxid=7 ! Pinus strobus
iF(baum(i).EQ.24) ngroups(nlines)%taxid=10 ! Douglasie
IF (dm(i).eq.0) dm(i) = 0.5
IF (mhoe(i).eq.0) mhoe(i) = 1.0
IF (gf(i).eq.0) gf(i) = 0.25
......@@ -310,50 +310,6 @@ if(infile=='input/hyyti_ini_0616.txt') then
3333 CONTINUE
nlines=nlines-1
WRITE(*,*) nlines,'sets of data', numstand, 'sets of stands'
ELSE
IF(select_lines) THEN
READ(listunit,*)nlines_comp
ALLOCATE(locid_comp(nlines_comp))
DO i=1,nlines_comp ! reading list of sites to be initialised
READ(listunit,*) locid_comp(i)
ENDDO ! end reading list of sites to be initialised
ENDIF ! end of reading file with sites to be selected
IF(select_lines) CLOSE(listunit)
CALL assign_DSW
CALL init_plenter_param
READ (inunit,*)nlines
ALLOCATE(ngroups(nlines))
istart=1
READ(inunit,*) ngroups(1)%locid,ngroups(1)%schicht,ngroups(1)%BRAid,ngroups(1)%alter,ngroups(1)%patchsize,ngroups(1)%mhoe,ngroups(1)%dm,ngroups(1)%volume,ngroups(1)%gf
ngroups(1)%patchsize=ngroups(1)%patchsize*10000.
ngroups(1)%baumzahl=0
ngroups(istart)%standsize=ngroups(1)%patchsize
ngroups(1)%taxid=tax_of_BRA_id(ngroups(1)%BRAid)
DO i=2,nlines
READ(inunit,*) ngroups(i)%locid,ngroups(i)%schicht,ngroups(i)%BRAid,ngroups(i)%alter,ngroups(i)%patchsize,ngroups(i)%mhoe,ngroups(i)%dm,ngroups(i)%volume,ngroups(i)%gf
WRITE(*,*) 'set no', i, 'read'
ngroups(i)%baumzahl=0
! the following line maps BRAid 770 to 779, other 'Mehlbeeren', because two
! different numbering systems existed in Brandenburg in the course of time
IF(ngroups(i)%BRAid==770) ngroups(i)%BRAid=779
ngroups(i)%patchsize=ngroups(i)%patchsize*10000.
ngroups(i)%taxid=tax_of_BRA_id(ngroups(i)%BRAid)
IF(ngroups(i)%taxid==6) ngroups(i)%taxid=3
IF(ngroups(i)%taxid==0) THEN
ELSE
ENDIF
IF(ngroups(i)%locid==ngroups(istart)%locid) THEN
ngroups(istart)%standsize=ngroups(istart)%standsize+ngroups(i)%patchsize
ngroups(i)%standsize = ngroups(istart)%standsize
ELSE
istart=i
ngroups(istart)%standsize=ngroups(i)%patchsize
fl_num=fl_num+1
ENDIF
ENDDO ! readin loop for multi data-set
ENDIF ! block for direct DSW data or brb_inv-file structure
CLOSE(inunit)
! read in file headder for description, write into ini-file
cform=1;hlp_lai=0
......@@ -571,6 +527,7 @@ if(infile=='input/hyyti_ini_0616.txt') then
! classification of single values in diameter cohorts
clwdth=1+AINT((bhdmax-bhdmin)/ncl1) !calculation of class widths
! write(4444,*) 'clwdth', clwdth, bhdmax, bhdmin, ncl1
DO i=1,ncl1
nz(i)=0
zbhd(i)=0
......@@ -771,8 +728,9 @@ CASE(6)
g=ngroups(iz)%gf !basal area/ha
gpatch=g*4. !basal area/patch
bz=ngroups(iz)%baumzahl*4. !tree numbre/patch
clwdth=dg/20.
! clwdth=dg/20.
clwdth=dg/5
! selection of uni-height curve: beech, spruce, oak calculation according to Weimann,
! other species of trees after Kuleschis (vergl. Gerold 1990)
IF (taxid==3.OR.taxid==5) THEN
......
......@@ -165,6 +165,7 @@ do
case (8, 9, 10)
call readsoil ! reading soil parameter
IF (flag_end .gt.0) return
call readredN ! Input redN or test resp.
end select
endif
......@@ -198,7 +199,7 @@ call readlit
! Initialization of soil model with profile data
call soil_ini ! Aufruf ohne s_cn_ini
! Initialization disturbances
IF (flag_dis .eq. 1) CALL dist_ini
IF (flag_dis .eq. 1 .or. flag_dis .eq. 2) CALL dist_ini
! Initialization of stand
call prepare_stand
IF (flag_end .gt.0) return
......@@ -268,12 +269,13 @@ if (flag_eva .gt.10) call evapo_ini
subroutine readsoil ! Input of soil parameter
use data_par
use data_soil_t
use data_site
implicit none
integer :: inunit, helpnl, helpnr
integer :: inunit, helpnl, helpnr, ihelp
real helpgrw, hlong, hlat
character :: text
character(30) :: hor, boart, helpid
......@@ -283,17 +285,15 @@ if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' readsoil'
! Setting of flag_surf from flag_cond
select case (flag_cond)
case (0,1,2,3)
flag_surf = 0
case (10,11,12,13)
flag_surf = 1
case (0,1,2,3)
flag_surf = 0
case (20,21,22,23)
flag_surf = 2
case (10,11,12,13)
flag_surf = 1
case (30,31,32,33)
flag_surf = 3
case (30,31,32,33)
flag_surf = 3
end select
! Setting of flag_bc from flag_decomp
......@@ -479,7 +479,7 @@ IF (ex .eqv. .true.) then
if (.not.flag_mult8910) print *,' >>>FORESEE message: soil_id ', soilid(ip), ' not found'
if (.not.flag_mult8910) print *,' Check your input choice!!!'
if (help==1) call dealloc_soil
CALL error_mess(time,"soil identificator not found"//adjustl(soilid(ip))//"ip No.",real(help_ip))
CALL error_mess(time,"soil identificator not found "//adjustl(soilid(ip))//"ip No.",real(help_ip))
flag_end = 5
return
ENDIF ! ios
......@@ -535,15 +535,11 @@ IF (ex .eqv. .true.) then
endif
end do
IF (ios .ne.0) then
if (.not.flag_mult8910) print *,' >>>FORESEE message: Error during reading soil data!'
WRITE(*,'(A)',advance='no') ' Stop program (y/n)? '
read *, a
IF ( a .eq. 'y' .or. a .eq. 'Y') then
print *, ' STOP program!'
stop
endif
print *,' >>>FORESEE message: Error during reading soil data!'
print *, ' Program stopped!'
IF (help==1) call dealloc_soil
if (.not.flag_mult8910) print *,' Check your input choice!!!'
flag_end = 7
return
endif ! ios
exit
endif
......@@ -556,7 +552,7 @@ IF (ex .eqv. .true.) then
print *,' Check your input choice!!!'
endif
if (help==1) call dealloc_soil
CALL error_mess(time,"soil identificator not found"//adjustl(soilid(ip))//"ip No.",real(help_ip))
CALL error_mess(time,"soil identificator not found "//adjustl(soilid(ip))//"ip No.",real(help_ip))
flag_end = 5
return
ENDIF ! ios
......@@ -915,7 +911,7 @@ if (.not.flag_mult8910 .or. (flag_mult8910 .and. anh .eq. "1") .or. (flag_mult89
WRITE(unit_ctr,'(A66,I4)') 'Time step for photosynthesis calculations (days) - ns_pro: ',ns_pro
WRITE(unit_ctr,'(A66,I4)') 'Mortality (0-OFF,1-ON stress, 2- ON stress+intr) - flag_mort: ',flag_mort
WRITE(unit_ctr,'(A66,I4)') 'Regeneration (0-OFF,1-ON, 2-weekly growth of seedl.) - flag_reg: ',flag_reg
WRITE(unit_ctr,'(A66,I4)') 'use FORSKA for regeneration (0-OFF,1-ON) - flag_forska: ',flag_forska
WRITE(unit_ctr,'(A66,I4)') 'use FORSKA for regeneration (0-OFF,1-ON) - flag_forska: ',flag_lambda
WRITE(unit_ctr,'(A66,I4)') 'Stand initialization (0-no,1-from *.ini,2-generate) - flag_stand: ',flag_stand
WRITE(unit_ctr,'(A66,I4)') 'Ground vegetation initialization (0-no,1-generate) - flag_sveg: ',flag_sveg
WRITE(unit_ctr,'(A66,I4)') 'Stand management (0-no,1-yes, 2 - seed once) - flag_mg: ',flag_mg
......@@ -1179,22 +1175,17 @@ real hNO, hNH
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' readdepo'
if (.not.allocated(NOd)) allocate (NOd (1:366,1:year))
if (.not.allocated(NHd)) allocate (NHd (1:366, 1:year))
! for areal usage standard/constant deposition is set as concentration:
if (flag_multi .eq. 8 .or. flag_mult910) then
flag_depo = 2
if (.not.allocated(NOd)) then
allocate (NOd (1:366,1:year))
NOd = NOdep(ip) ! concentration mg/l
endif
if (.not.allocated(NHd)) then
allocate (NHd (1:366,1:year))
NHd = NHdep(ip) ! concentration mg/l
endif
NOd = NOdep(ip) ! concentration mg/l
NHd = NHdep(ip) ! concentration mg/l
return
endif
if (.not.allocated(NOd)) allocate (NOd (1:366,1:year))
if (.not.allocated(NHd)) allocate (NHd (1:366, 1:year))
NOd = 0.
NHd = 0.
......@@ -1427,6 +1418,7 @@ END subroutine readdepo
SUBROUTINE readredN
use data_out
use data_site
use data_species
use data_stand
use data_simul
......
......@@ -32,7 +32,6 @@ character a
character(8) actdate
character(10) acttime, helpsim, text1, text2
real time1, time2, time3
logical lhelp
unit_err=getunit()
if(flag_multi.eq.5) dirout = './'
......@@ -46,20 +45,6 @@ write (unit_trace, '(I4,I10,A)') iday, time_cur, ' sim_control'
! check daily output
if (year > 5 .and. flag_dayout .ge. 1) then
lhelp = .true.
do i = 1,outd_n
if (outd(i)%out_flag .eq. flag_dayout) then
select CASE (outd(i)%kind_name)
CASE ('day_short')
lhelp = .false.
end select
endif
enddo
if (lhelp) 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)? '
......@@ -67,7 +52,6 @@ if (year > 5 .and. flag_dayout .ge. 1) then
IF (a .eq. 'n' .or. a .eq. 'N') then
flag_dayout = 0
ENDIF
endif ! lhelp
ENDIF
! open file ycomp (yearly compressed output (multi run))
......@@ -136,10 +120,14 @@ time3 = 0.
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
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
case default
print*,ip, 'flag_end = ', flag_end
end select
......@@ -228,7 +216,7 @@ if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' simulation_4C
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) then
if (flag_standup .gt. 0 .or. flag_dis==1 .or. flag_dis==1) then
call stand_balance
call standup
flag_standup = 0
......@@ -239,7 +227,10 @@ if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' simulation_4C
! read or create Redn for areal application
IF (time.EQ.1 .and. flag_redn) CALL RedN_ini
IF (flag_dis .eq. 1) CALL dis_manag
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
......
......@@ -453,12 +453,6 @@ CASE (0, 10, 20, 30, 40) ! de Vries
tcond0 = numera/denom * 86400. ! s --> day
CASE(2, 12, 22, 32, 42) ! sum like resistor; wie Widerstaende addieren
tcond2 = water%vf / water%tc + quarz%vf / quarz%tc + clay%vf / clay%tc + &
silt%vf / silt%tc + humus%vf / humus%tc + air%vf / air%tc + stone%vf / stone%tc + ice%vf / ice%tc
tcond2 = 86400. / tcond2
CASE(3, 13, 23, 33, 43) ! Campbell
vfm = clay%vf + silt%vf + stone%vf
vfs = vfm + quarz%vf + humus%vf
......@@ -504,13 +498,6 @@ CASE (1, 11, 21, 31, 41) ! Neusypina
hcapi = hcap1
tcondi = tcond1
CASE (2, 12, 22, 32, 42) ! sum like resitors; Widerstnde addieren
if ((tcond2 .gt. 8000.) .or. (tcond2 .le. 0.)) then
continue
endif
hcapi = hcap0
tcondi = tcond2
CASE (3, 13, 23, 33, 43) ! Campbell
hcapi = hcap0
tcondi = tcond3
......
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Subroutines used only with flag flag_forska
!
! cetbl_4c
! CGTSPE_4c
! CLIMEF_4c
! gsdr_cal
! tmp_mean
! therm
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE CETBL_4c
use data_effect
use data_taxa
use data_simul
use data_stand
! function declarations
REAL RAND
! local variables
real :: PMX
INTEGER :: I,J,K
integer,dimension(17) :: nsap= 0
real,dimension(17) :: amdest = 0., &
amdest1 = 0.
if (flag_light.eq.1.or.flag_light.eq.2) then
PMX= Vstruct(lowest_layer)%Irel
else if (flag_light.eq.3.OR.flag_light.EQ.4) then
PMX = Bgpool(lowest_layer+1)
end if
! amend the EST for climate according to the climate multipliers
do i=1,17
AMDEST(I)=EST(I)*GDDMX(I)*DRMX(I)*TCMX(I)*TWMX(I)*PMX &
*XTFTMX(I)*TWARMX(I)
AMDEST1(I)=EST(I)*AMIN1(GDDMX(I),DRMX(I),TCMX(I),TWMX(I), &
PMX,XTFTMX(I),TWARMX(I))
IF(GSC(I).EQ.0.0)GOTO 301
301 CONTINUE
end do
RETURN
END subroutine cetbl_4c
SUBROUTINE CGTSPE_4c
! input of species data for regeneration
! reads species parameters
use data_simul
use data_taxa
! local variables
INTEGER:: I,J,K,nowunit,ntax
! reads number of taxa (NTAX)
nowunit=getunit()
open(unit=nowunit,file= '/data/safe/4C/4C_input/par/param_4c.dat', status='old')
READ(nowunit,*) NTAX
! reads for each taxon:
! NAM(I): name (8 characters)
! HMX(I): max height (m)
! HDS(I): initial slope of diameter vs height (m/cm)
! hgro(I): maximum height growth per year (m)
! ALP(I): half-saturation point (umol/m**2/s)
! LCP(I): compensation point (umol/m**2/s)
! GSC(I): growth constant (cm**2/m/yr)
! EST(I): sapling establishment rate (/ha/yr)
! TDI(I): threshold relative growth efficiency for increased mortality
! UMN(I): intrinsic mortality rate (/yr)
! UMX(I): suppressed mortality rate (/yr)
! SPR(I): number of sprouts per tree (0.0 or greater)
! SMN(I): minimum diameter for sprouting (cm)
! LAC(I): initial leaf area/D2 ratio (m**2/cm**2)
! LAF(I): sapwood turnover rate (/yr)
! BCF(I): stemwood biomass conversion factor (kg/cm**2/m)
! R(I): volumetric sapwood maintenance cost (/yr)
! Q10(I): rate of increase of respiration
! TMIN(I): minimum temperature for assimilation
! TMAX(I): maximum temperature for assimilation
! CCP(I): species compensation point
! DRI(I): maximum tolerated drought-index
!MINGDD(I): minimum growing degree-days
! MINTC(I): minimum temperature of coldest month (degrees C)
! MAXTC(I): maximum temperature of coldest month (degrees C)
! MINTW(I): minimum temperature of warmest month (degrees C)
! DORE(I): deciduous or evergreen 0=deciduous,1=evergreen
! ntc(I): nitrogen tolerance class (1,2,3,4,5)
! e1(I): Parameter smin of haadee height growth function
! e2(I): Second Parameter of haadee height growth function
! geff(I): growth efficiency factor of shaded trees
DO I=1,ntax
READ(nowunit,1) NAM(I)
READ(nowunit,*) HMX(I),HDS(I),hgro(I),ALP(I),LCP(I),GSC(I), &
EST(I),TDI(I),UMN(I),UMX(I),SPR(I),SMN(I),LAC(I),LAF(I),BCF(I), &
R(I),Q10(I),TMIN(I),TMAX(I),CCP(I),DRI(I),MINGDD(I),MINTC(I), &
MAXTC(I),MINTW(I),DORE(I),ntc(I)
IF(SPR(I).EQ.0)SMN(I)=0.0
DRI(I)=DRI(I)+0.3
end do
RETURN
! format statements
1 FORMAT(A8)
END subroutine cgtspe_4c
SUBROUTINE CLIMEF_4c
use data_taxa
use data_effect
use data_simul
! computes the growth multipliers.
! checks to see if GDD, temp coldest month below minimum for species
! if so multipliers = 0 else equals 1.
! computes drought effect multipliers as per ICP
! sets max.temp of coldest month multiplier to 0 or 1 for ESTBL routine
! checks if warmest month exceeds species limit
! averages light intensity (INS) over time step.
! local parameters
INTEGER :: I,J,K
REAL ::TOTGDD= 0, &
TGSDRT=0., &
TM4DRT=0.
real,dimension(17) :: tottft=0.
! gives growth multiplier for each species to be applied in subroutine
! TVXT or ETBL - growing degree days, growing/-4 drought index, temps.
TOTGDD=GDD(time)
TGSDRT=GSDRI(time)
TM4DRT=M4DRI(time)
! totals and then averages species specific multipliers etc. over timestep
! that is sapres, mutmx, tftmx
do i=1,17
xtftmx(i) = tftmx(i,time)
end do
! set multipliers to 1 before checking on environment
do i=1,17
GDDMX(I)=1.0
TWARMX(I)=1.0
TCMX(I)=1.0
TWMX(I)=1.0
TWARMX(I)=1.0
! check to see is a deciduous species
IF(DORE(I).EQ.0)THEN
DRMX(I)=1-((TGSDRT/DRI(I))**2)
IF(DRMX(I).LT.0.0)DRMX(I)=0.0
ELSE
! must be an evergreen
DRMX(I)=1-((TM4DRT/DRI(I))**2)
IF(DRMX(I).LT.0.0)DRMX(I)=0.0
ENDIF
! check if environment exceeds species limits - step functions
! if so set multiplier to zero
IF(TOTGDD.LT.MINGDD(I))GDDMX(I)=0.0
IF(TCOLD.LT.MINTC(I))TCMX(I)=0.0
IF(TCOLD.GT.MAXTC(I))TWMX(I)=0.0
IF(TWARM.LT.MINTW(I))TWARMX(I)=0.0
! write out to screen and forcli.out multipliers for each species
! keep these commented as they use a lot of paper <--M.B was ist damit gemeint? ist das relevant fr den nutzer.
end do
do i=1,17
end do
end subroutine climef_4c
SUBROUTINE gsdr_cal
! calculation of gsdri and m4dri for FORSKA regeneration
use data_climate
use data_effect
use data_simul
use data_evapo
if(tp(iday,time).ge.-4.) then
foudpt = foudpt + pet
foudae = foudae + aet
end if
if(tp(iday,time).ge.4.) then
tgsdpt = tgsdpt + pet
tgsdae = tgsdae + aet
end if
if(iday.eq. recs(time)) then
gsdri(time) = (tgsdpt-tgsdae)/tgsdpt
m4dri(time) = (foudpt-foudae)/foudpt
end if
END SUBROUTINE gsdr_cal
SUBROUTINE tmp_mean
! calculation of environmental variables twarm, tcold and long-term monthly
! mean of temperature
USE data_effect
USE data_climate
USE data_simul
real,dimension(12) :: tmph = 0.
integer :: i,l,m,dayc
allocate( tpmean(12))
allocate (gdd(year))
allocate (tftmx(17,year))
monrec=(/31,28,31,30,31,30,31,31,30,31,30,31/)
tpmean = 0
if (recs(time).eq.366) then
monrec(2)=29
else
monrec(2)=28
endif
do k = 1, year
! call calculation of env. variables
call therm(k)
dayc = 1
do l= 1,12
tmph(l) = 0.
do m=1,monrec(l)
tmph(l) = tmph(l) + tp( dayc,k)
dayc = dayc + 1
end do
tmph(l) = tmph(l)/monrec(l)
tpmean(l) = tpmean(l) + tmph(l)
end do
end do
do l=1,12
tpmean(l) = tpmean(l)/year
end do
! work out which is temperature of coldest month
! and warmest month for year
tcold = 50.0
twarm = -50.0
do k=1,12
if(tpmean(k).lt.tcold) tcold = tpmean(k)
if(tpmean(k).gt.twarm) twarm = tpmean(k)
end do
END SUBROUTINE tmp_mean
SUBROUTINE therm(ktime)
! therm - calculation of environmental variables (annual and species specific)
! gdd - growing degress day
! tftmx - thermal multiplier - species specific
use data_climate
use data_simul
use data_effect
use data_taxa
implicit none
! local variables
integer :: j,k,m4day,gdday1,ktime
real,dimension(17) :: tft,tresft
gdd(ktime) = 0.
m4day=0
gdday1=0
do j=1,17
tft(j)=0.0
tresft(j)=0.0
end do
! calculate ft values for each day of the year
! for each species upto number of taxa
do k=1,17
do j=1,recs(ktime)
! add up mutmx multiplier
tresft(k) = tresft(k)+(q10(k)**((tp(j,ktime) - tref)*0.1))
if(k.eq.1) then
if (tp(j,ktime).ge.tref) gdd(ktime) = gdd(ktime) + (tp(j,ktime)-tref)
end if
! first check to see if deciduous or not
if(dore(k).eq.0)then
! totalling daily deciduous multipliers for growing season only
if(tp(j,ktime).ge.5.0) then
tft(k) = tft(k)+(4*(tp(j,ktime)-tmin(k))*(tmax(k)-tp(j,ktime))/(tmin(k)-tmax(k))**2)
endif
else
! must be evergreen so produce daily values
! do not allow below zero
! checks for temperature greater than -4 oC for evergreen species
if(tp(j,ktime).ge.-4.0)then
tft(k)=tft(k)+(4*((tp(j,ktime)-tmin(k))*(tmax(k)-tp(j,ktime))) &
/(tmin(k)-tmax(k))**2)
endif
endif
if(tft(k).lt.0.0)tft(k)=0.0
end do
end do
do j=1,recs(ktime)
if(tp(j,ktime).ge.5.0) then
gdday1=gdday1+1
end if
if(tp(j,ktime).ge.-4.0) then
m4day=m4day+1
end if
end do
do k=1,17
if(dore(k).eq.0) then
tftmx(k,ktime) = tft(k)/gdday1
else
tftmx(k,ktime) = tft(k)/m4day
end if
end do
END SUBROUTINE therm
......@@ -104,7 +104,12 @@ do
sumbio = sumbio + ntr * zeig%coh%totBio
sumNPP = sumNPP + ntr * zeig%coh%NPP
Ndem = Ndem + ntr * zeig%coh%Ndemc_c
autresp = autresp + ntr * zeig%coh%maintres
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
......@@ -498,10 +503,6 @@ svar%frt = 0.
END SELECT
Enddo
IF(spar(spec_new)%phmodel==4) THEN
svar(spec_new)%daybb = svar(spec_new)%ext_daybb
ELSE
......@@ -543,7 +544,7 @@ do i=1,nspecies
svar(i)%mean_diam = svar(i)%mean_diam / ntr
svar(i)%mean_jrb = svar(i)%mean_jrb / ntr
svar(i)%basal_area = pi*ntr*(svar(i)%med_diam*svar(i)%med_diam/40000)*10000/kpatchsize
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
......@@ -1083,8 +1084,14 @@ IMPLICIT NONE
DO
IF(.not.associated(zeig)) exit
if (zeig%coh%species.ne.nspec_tree+2) then ! exclude mistletoe from senescence
zeig%coh%sfol = spar(zeig%coh%species)%psf * zeig%coh%x_fol
zeig%coh%sfrt = spar(zeig%coh%species)%psr * zeig%coh%x_frt
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
......
SUBROUTINE win_prepini
!*
!* Implementation and Revisions:
!* -----------------------------
!*
! 15.08.11 Su Zwischendruck bei treeini raus
! 24.02.11 Su stand_id als Character
! 01.02.05 ag win_prepini: help_ip set to actual run number
! 28.01.05 Su
! check stand_id from multi-stand-intialisaton file (...tree.ini)
! lmulti = 0 ==> single initialisation file
USE data_simul
USE data_site
USE data_stand
USE data_species
use flag_field
IMPLICIT NONE
CHARACTER :: a
CHARACTER(30) :: filename, text
CHARACTER(50) :: test_stand_id
INTEGER :: ios,treeunit
LOGICAL :: exs, lstin
INTEGER :: help_ip, test_vf, i
REAL :: test_patchsize, xx
IF(site_nr==1) THEN
help_ip=site_nr
ELSE
help_ip=act_run
END IF
! when initialization stand flag == 1
IF(flag_stand>0) then
exs = .false.
stand_id = standid(help_ip)
! reading stand information from treefile
inquire (File = treefile(help_ip), exist = exs)
IF(exs .eqv. .false.) write(*,*) ' Stand initialization file not exists!'
! read values from treefile
IF (exs.eqv. .true.) then
treeunit=getunit()
OPEN(treeunit,file=treefile(help_ip),action='read', pad='YES')
READ(treeunit,'(I1,F12.0)',iostat=ios) test_vf, test_patchsize
IF(test_patchsize .GT. 0.) THEN
lmulti = .FALSE.
IF(test_patchsize.NE.kpatchsize) THEN
CALL error_mess(time,"patch size in sim-file and the one used for initialisation do not match",kpatchsize)
CALL error_mess(time,"value in ini-file",test_patchsize)
CALL error_mess(time,"value in sim-file",kpatchsize)
kpatchsize = test_patchsize
ENDIF
ELSE
lmulti = .TRUE.
! count number of stand_id
anz_standid = 0
do
READ(treeunit,'(A)',iostat=ios) a
IF (a .ne. '!') exit
end do
do
read(treeunit,*,iostat=ios) xx
IF (ios .lt. 0) exit
anz_standid = anz_standid + 1
do while (xx .gt. -90.0)
read (treeunit,*) xx
! print *,xx
enddo ! xx
enddo
if (allocated(standid_list)) deallocate(standid_list)
allocate (standid_list(anz_standid))
rewind treeunit
! read stand_id
read (treeunit,*) xx
do
READ(treeunit,'(A)',iostat=ios) a
IF (a .ne. '!') exit
end do
backspace treeunit
lstandid = .FALSE.
do i = 1, anz_standid
read(treeunit,*) test_stand_id, test_patchsize, text
standid_list(i) = test_stand_id
if (test_stand_id .eq. stand_id) lstandid = .TRUE.
read (treeunit,*) xx
do while (xx .gt. -90.0)
read (treeunit,*) xx
enddo ! xx
enddo
END IF ! test_patchsize -- lmulti
CLOSE(treeunit)
END IF ! exs
END IF
END SUBROUTINE win_prepini
......@@ -154,7 +154,7 @@ integer act_spec, act_year, set_year
! calculate emission from harvesting process
emission_har (management_years(i)) = summe * sub_par(1)
write (9999,*) emission_har(management_years(i)), management_years(i)
! write (9999,*) emission_har(management_years(i)), management_years(i)
end do
end subroutine calculate_product_lines
......