!*****************************************************************!
!* *!
!* 4C (FORESEE) *!
!* *!
!* *!
!* Subroutines: *!
!* PREPARE_SITE and PREPARE_CLIMATE *!
!* *!
!* Contains subroutines: *!
!* *!
!* PREPARE_SITE: *!
!* preparation of site specific simulation parameters *!
!* *!
!* contains internal subroutines: *!
!* SITEMENU: choice of inputs *!
!* EDITFILE: edit filenames *!
!* READSOIL: input of soil parameter *!
!* READCN: input of C-N-parameter *!
!* READVALUE: input of start values for *!
!* soil water and C-N-modeling *!
!* ALLOC_SOIL: allocate soil variables *!
!* STAND_BAL_INI: allocate and init stand variables *!
!* CONTROL_FILE: saving all parameters *!
!* and start conditions for each site *!
!* *!
!* READDEPO: reading deposition data *!
!* READREDN: reading values of redN *!
!* READLIT: reading initialisation data of litter fractions *!
!* *!
!* PREPARE_CLIMATE: reading of site specific climate input data *!
!* from file *!
!* contains internal subroutines: *!
!* READ_DWD *!
!* READ_CLI *!
!* CLIMFILL *!
!* *!
!* STORE_PARA: multi run - restore of changed parameter *!
!* *!
!* 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 prepare_site
! input of site specific data
use data_climate
use data_inter
use data_manag
use data_mess
use data_out
use data_par
use data_simul
use data_site
use data_soil
use data_soil_cn
use data_species
use data_stand
use data_tsort
use data_frost
implicit none
integer i,ios,help, help_ip
character a
character :: text
character(10) :: helpsim, text2
logical:: ex=.TRUE.
real parerr
real, external :: avg_sun_incl
character(100) :: helpx
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' prepare_site'
WRITE(helpsim,'(I4)') anz_sim
read(helpsim,*) anh
IF(site_nr==1) THEN
help_ip=site_nr
ELSE
help_ip=ip
END IF
! Initialization of climate data
IF (flag_clim==1 .or. ip==1 .or. flag_multi .eq.5) THEN
call prepare_climate
END IF
if (flag_end .gt. 0) return
ios=0; help=0
do
if (ip==1 .and. flag_mult9) then
if (flag_trace) write (unit_trace, '(I4,I10,A,I3,A5,L5)') iday, time_cur, ' prepare_site ip=',ip,' ex=',ex
call readspec
call readsoil ! reading soil parameter
IF (flag_end .gt.0) return
if (flag_soilin .eq. 0) call readvalue ! Initialization of simulation start values for soil layers
! biochar
if (flag_bc .gt. 0) call bc_appl
! Deposition data
call readdepo
! Input redN
if (flag_multi .ne. 4 .or. flag_multi .ne. 8 ) call readredN
flag_mult9 = .FALSE.
else
if (flag_trace) write (unit_trace, '(I4,I10,A,I3,A5,L5)') iday, time_cur, ' prepare_site ip=',ip,' ex=',ex
! Deposition data
call readdepo
select case (flag_multi)
case (1,6)
call readspec
if (flag_soilin .eq. 0) call readvalue ! Initialization of simulation start values for soil layers
call readredN ! Input redN
call readsoil ! reading soil parameter
do
jpar = jpar + 1
if (vpar(jpar) .gt. -90.0) then
helpx = simpar(jpar)
call store_para(vpar(jpar), helpx, parerr)
IF (parerr .ne. 1.) then
CALL error_mess(time,'parameter variation: '//trim(simpar(jpar))//' not found',vpar(jpar))
write (*,*) '*** parameter variation: ', trim(simpar(jpar)), ' not found, see error log'
endif
else
exit
endif
enddo
! biochar
if (flag_bc .gt. 0 .or. flag_decomp .gt. 100) call bc_appl
case (2,4)
call readsoil ! reading soil parameter
if (flag_soilin .eq. 0) call readvalue ! Initialization of simulation start values for soil layers
case (5)
call readspec
call readsoil
if (flag_soilin .eq. 0) call readvalue ! Initialization of simulation start values for soil layers
call readredN ! Input redN
case (7)
call assign_co2par
call readsoil ! reading soil parameter
if (flag_soilin .eq. 0) call readvalue ! Initialization of simulation start values for soil layers
call readredN ! Input redN
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
exit
enddo
! Setting flag_inth and prec_stad_red from flag_int
if (flag_int .lt. 1000) then
flag_inth = flag_int
else
! Conversion character ==> number and vice versa
write (helpsim,'(I4)') flag_int
text2 = helpsim(2:2)
read (text2,*) flag_inth
text2 = helpsim(3:4)
read (text2,*) prec_stand_red
endif
if (.not.flag_mult8910) then
unit_soil = getunit()
open (unit_soil,file=trim(dirout)//trim(site_name(help_ip))//'_soil.ini'//anh,status='replace')
WRITE (unit_soil,'(2A)') '! Soil initialisation, site name: ',site_name(help_ip)
endif
call stand_bal_ini !allocation of stand summation variables
! Initialization of CO2
call assign_co2par
! Initialisation litter compartments
call readlit
! Initialization of soil model with profile data
call soil_ini ! Aufruf ohne s_cn_ini
! Initialization disturbances
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
! calculation of latitude in radians
xlat = lat/90.*pi*0.5
! calculation of average sun inclination
avg_incl = AVG_SUN_INCL(lat) ! degrees
beta=avg_incl*PI/180 ! radians
! read externally prescribed bud burst days
CALL readbudb
! Initialization management
IF(flag_mg.ne.0.and. flag_mg.ne.5) call manag_ini
IF(flag_mg.eq. 5) then
thin_dead = 1
allocate(thin_flag1(nspec_tree))
thin_flag1 = 0
end if
! Initialization of output file per site
call prep_out
call stand_balance
call CROWN_PROJ
call standup
call root_ini ! initialisation of root distribution
call s_cn_ini
! Initialization of soil temperature model with stand data
call s_t_ini
! control file for saving simulation environment
! output of first Litter-Input at start
if(flag_mult8910 .and. (anz_sim .gt. 1)) then
continue
else
IF ((ip .eq. 1 .or. flag_multi .eq. 1 .or. flag_multi .eq. 6) .and. (time_out .ne. -2) ) call control_file
endif
! hand over of the litter-initialising
call litter
if ((flag_decomp .eq. 20) .or. (flag_decomp .eq. 21)) then
call testfile(valfile(ip),ex)
if (ex .eqv. .true.) then
ios = 0
unit_litter = getunit()
open(unit_litter,file=valfile(ip),status='old',action='read')
if (flag_multi .ne. 9) print *,' *** Open file of litter input data ',valfile(ip),'...'
do
read(unit_litter,*) text
IF(text .ne. '!')then
backspace(unit_litter);exit
endif
enddo
endif
endif
call cn_inp
! read flux data
if (flag_eva .gt.10) call evapo_ini
! yearly output
IF (time_out .gt. 0) THEN
IF (mod(time,time_out) .eq. 0) CALL outyear (1)
IF (mod(time,time_out) .eq. 0) CALL outyear (2)
ENDIF
contains
!-------------------------------------------------------------------------------
subroutine readsoil ! Input of soil parameter
use data_par
use data_soil_t
use data_site
implicit none
integer :: inunit, helpnl, helpnr, ihelp
real helpgrw, hlong, hlat
character :: text
character(30) :: hor, boart, helpid
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 (30,31,32,33)
flag_surf = 3
end select
! Setting of flag_bc from flag_decomp
if (flag_decomp .ge. 100) then
flag_decomp = flag_decomp - 100
flag_bc = 1
else
flag_bc = 0
endif
call testfile(sitefile(ip),ex)
IF (ex .eqv. .true.) then
inunit = getunit()
ios=0
open(inunit,file=sitefile(ip),iostat=ios,status='old',action='read')
if (.not.flag_mult8910) then
print *,'***** Reading soil parameter from file ',sitefile(ip),'...'
write (unit_err, *) 'Soil parameter from file ',trim(sitefile(ip))
endif
do
read(inunit,*) text
IF(text .ne. '!')then
backspace(inunit)
exit
endif
enddo
if (flag_multi .eq. 8.or. flag_multi.eq.5.or. flag_mult910) then
read(inunit,*) text
IF((text .eq. 'N') .or. (text .eq. 'n'))then
flag_soilin = 3
else
flag_soilin = 2
backspace(inunit)
endif
else
read(inunit,*) text
IF((text .eq. 'N') .or. (text .eq. 'n'))then
flag_soilin = 1
else
flag_soilin = 0
backspace(inunit)
endif
soilid(ip) = valfile(ip)
endif
if ((text .eq. 'S') .or. (text .eq. 's'))then
flag_soilin = 4
read(inunit,*) text
endif
if (.not.flag_mult8910) then
write (unit_err, *) 'Soil identity number ', trim(soilid(ip))
write (unit_err, *) 'Climate ID ', trim(clim_id(ip))
endif
if (flag_soilin .eq. 1 .or. flag_soilin .ge. 3) then
flag_hum = 1
endif
if (flag_cond .ge. 40) then
flag_hum = 0
endif
select case (flag_soilin)
case (0,1) ! single files f. j. site
read (inunit,*,iostat=ios) long
read (inunit,*,iostat=ios) lat
read (inunit,*,iostat=ios) nlay
read (inunit,*,iostat=ios) nroot_max
read (inunit,*,iostat=ios) helpgrw
if (helpgrw .gt. 1) then
grwlev = helpgrw
else
fakt = helpgrw
grwlev = 1000.
endif
read (inunit,*,iostat=ios) w_ev_d
read(inunit,*,iostat=ios) k_hum ! mineralization constants of humus
read(inunit,*,iostat=ios) k_hum_r
read(inunit,*,iostat=ios) k_nit ! nitrification constant
IF(help==0) call alloc_soil
read (inunit,*,iostat=ios) text
select case (flag_soilin)
case (0) ! old input structure
do i = 1, nlay
read (inunit,*,iostat=ios) text
read (inunit,*,iostat=ios) thick(i),pv_v(i),dens(i),f_cap_v(i), &
wilt_p_v(i),spheat(i),phv(i),wlam(i)
end do
skelv = 0.
case(1) ! new input structure
do i = 1, nlay
read (inunit,*,iostat=ios) helpnr, thick(i),pv_v(i),f_cap_v(i),wilt_p_v(i), &
dens(i),spheat(i),phv(i),wlam(i),skelv(i), sandv(i),clayv(i),humusv(i),&
C_hum(i), N_hum(i),NH4(i),NO3(i)
if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
if (phv(i) .le. 0.01) phv(i)=6.0 ! if flag_wurz 4 or 6 is used for calculation a pH-value is assumed
endif
end do
end select ! flag_soilin (0,1)
if (.not.flag_mult8910) print *, ' '
IF (ios .ne.0) then
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
IF (help==1) call dealloc_soil
print *,' Check your input choice!!!'
endif ! ios
case (2) ! all sites are read from one file; old structure
ios = 0
do while (ios .eq. 0)
read (inunit,*,iostat=ios) helpid, helpnl, helpnr
if (trim(soilid(ip)) .ne. trim(helpid)) then
do i = 1, helpnl
read (inunit,*,iostat=ios) helpid
enddo
else
nlay = helpnl
nroot_max = helpnr
if (help==0) call alloc_soil
do i = 1, nlay
read (inunit,*,iostat=ios) helpnl, hor, boart, depth(i), thick(i),pv_v(i),dens(i), &
f_cap_v(i), wilt_p_v(i), spheat(i),phv(i),wlam(i), &
C_hum(i), N_hum(i), NH4(i), NO3(i), temps(i)
enddo
lat = latitude(ip)
grwlev = gwtable(ip)
exit
endif
enddo
IF (ios .lt. 0) 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))
flag_end = 5
return
ENDIF ! ios
skelv = 0.
case (3) ! all sites are read from one file; new structure
ios = 0
do while (ios .eq. 0)
read (inunit,*,iostat=ios) helpid, helpnl, helpnr
if (trim(soilid(ip)) .ne. trim(helpid)) then
do i = 1, helpnl
read (inunit,*,iostat=ios) helpid
enddo
else
nlay = helpnl
nroot_max = helpnr
if (help==0) call alloc_soil
do i = 1, nlay
read (inunit,*,iostat=ios) helpnr, hor, boart, depth(i), thick(i),pv_v(i),f_cap_v(i), &
wilt_p_v(i),dens(i),spheat(i),phv(i),wlam(i),skelv(i), sandv(i), &
clayv(i),humusv(i),C_hum(i), N_hum(i),NH4(i),NO3(i)
if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
if (phv(i) .le. 0.01) phv(i)=6.0 ! if flag_wurz 4 or 6 is used for calculation a pH-value is assumed
endif
end do
lat = latitude(ip)
grwlev = gwtable(ip)
exit
endif
enddo
IF (ios .lt. 0) 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))
flag_end = 5
return
ENDIF ! ios
case (4) ! one file several sites
if (.not.flag_mult8910) print *,' Reading soil model parameter from soil type file... ', soilid(ip)
ios = 0
do while (ios .eq. 0)
read (inunit,*,iostat=ios) helpid
if (trim(soilid(ip)) .ne. trim(helpid)) then
read (inunit,*,iostat=ios) text
read (inunit,*,iostat=ios) text
read (inunit,*,iostat=ios) helpnl
do i = 1, helpnl+6
read (inunit,*,iostat=ios) boart
enddo
read (inunit,*,iostat=ios) boart
else
read (inunit,*,iostat=ios) hlong
read (inunit,*,iostat=ios) hlat
read (inunit,*,iostat=ios) nlay
read (inunit,*,iostat=ios) nroot_max
read (inunit,*,iostat=ios) helpgrw
if (flag_multi .eq. 8.or. flag_multi.eq.5.or. flag_mult910) then
if (abs(latitude(ip)) .gt. 90.) lat = latitude(ip)
grwlev = gwtable(ip)
else
if (helpgrw .gt. 1) then
grwlev = helpgrw
else
fakt = helpgrw
grwlev = 1000.
endif
long = hlong
lat = hlat
endif
read (inunit,*,iostat=ios) w_ev_d
read(inunit,*,iostat=ios) k_hum ! mineralization constants of humus
read(inunit,*,iostat=ios) k_hum_r
read(inunit,*,iostat=ios) k_nit ! nitrification constant
IF(help==0) call alloc_soil
read (inunit,*,iostat=ios) text
do i = 1, nlay
read (inunit,*,iostat=ios) helpnr, thick(i),pv_v(i),f_cap_v(i),wilt_p_v(i), &
dens(i),spheat(i),phv(i),wlam(i),skelv(i), sandv(i),clayv(i),humusv(i),&
C_hum(i), N_hum(i),NH4(i),NO3(i)
if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
if (phv(i) .le. 0.01) phv(i)=6.0 ! if flag_wurz 4 or 6 is used for calculation a pH-value is assumed
endif
end do
IF (ios .ne.0) then
print *,' >>>FORESEE message: Error during reading soil data!'
print *, ' Program stopped!'
IF (help==1) call dealloc_soil
flag_end = 7
return
endif ! ios
exit
endif
enddo
if (.not.flag_mult8910) print *, ' '
IF (ios .lt. 0) then
if (.not.flag_mult8910) then
print *,' >>>FORESEE message: soil_id ', soilid(ip), ' not found'
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))
flag_end = 5
return
ENDIF ! ios
end select ! flag_soilin
close(inunit)
endif ! ex
if (nroot_max .lt. 0) then
do i=1, nlay
if (C_hum(i) .gt. zero) nroot_max = i
enddo
endif
if (.not.flag_mult8910) then
write (unit_err, *) 'Latitude ',lat
write (unit_err,*)
endif
end subroutine readsoil
!-------------------------------------------------------------------------
subroutine readvalue ! Input of cn-parameters and start values for soil model
integer :: inunit
character :: text
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' readvalue'
call testfile(valfile(ip),ex)
IF (ex .eqv. .true.) then
ios = 0
inunit = getunit()
open(inunit,file=valfile(ip),status='old',action='read')
if (.not.flag_mult8910) print *,' *** Reading initial soil values from file ',valfile(ip),'...'
do
read(inunit,*) text
IF(text .ne. '!')then
backspace(inunit);exit
endif
enddo
! Soil temperature
read(inunit,*,iostat=ios) text
read(inunit,*,iostat=ios) (temps(i),i=1,nlay)
! C-content of humus
read(inunit,*,iostat=ios) text
read(inunit,*,iostat=ios) (C_hum(i),i=1,nlay)
! N-content of humus
read(inunit,*,iostat=ios) text
read(inunit,*,iostat=ios) (N_hum(i),i=1,nlay)
! NH4-content
read(inunit,*,iostat=ios) text
read(inunit,*,iostat=ios) (NH4(i),i=1,nlay)
! NO3-content
read(inunit,*,iostat=ios) text
read(inunit,*,iostat=ios) (NO3(i),i=1,nlay)
endif
IF (ios .ne. 0) then
print *,' >>>FORESEE message: Error during reading start values!'
WRITE(*,'(A)',advance='no') ' Stop program (y/n)? '
read *, a
IF ( a .eq. 'y' .or. a .eq. 'Y') then
print *, ' STOP program!'
stop
ELSE
call dealloc_soil
print *,' Check your input choice!!!'
end if
endif
close(inunit)
end subroutine readvalue
!--------------------------------------------------------------------------
subroutine alloc_soil
use data_soil_t
use data_soil
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' alloc_soil'
help=0
allocate(thick(nlay))
allocate(mid(nlay))
allocate(depth(nlay))
allocate(pv(nlay))
allocate(pv_v(nlay))
allocate(dens(nlay))
allocate(f_cap_v(nlay))
allocate(field_cap(nlay))
allocate(wilt_p(nlay))
allocate(wilt_p_v(nlay))
allocate(vol(nlay))
allocate(quarzv(nlay))
allocate(sandv(nlay))
allocate(BDopt(nlay))
allocate(clayv(nlay))
allocate(siltv(nlay))
allocate(humusv(nlay))
allocate(fcaph(nlay))
allocate(wiltph(nlay))
allocate(pvh(nlay))
allocate(dmass(nlay))
allocate(skelv(nlay))
allocate(skelfact(nlay))
allocate(spheat(nlay))
allocate(phv(nlay))
allocate(wlam(nlay))
allocate(wats(nlay))
allocate(watvol(nlay))
allocate(wat_res(nlay))
wat_res = 0.
allocate(perc(nlay))
allocate(wupt_r(nlay))
allocate(wupt_ev(nlay))
allocate(s_drought(nlay))
allocate(root_fr(nlay))
!allocate(dp_rfr(nlay))
allocate(temps(nlay))
allocate (C_opm(nlay))
allocate (C_hum(nlay))
allocate (C_opmfrt(nlay))
allocate (C_opmcrt(nlay))
allocate (N_opm(nlay))
allocate (N_hum(nlay))
allocate (N_opmfrt(nlay))
allocate (N_opmcrt(nlay))
allocate (NH4(nlay))
allocate (NO3(nlay))
allocate (Nupt(nlay))
allocate (Nmin(nlay))
allocate (rmin_phv(nlay))
allocate (rnit_phv(nlay))
allocate (cnv_opm(nlay))
allocate (cnv_hum(nlay))
allocate(slit(nspecies))
allocate(slit_1(nspecies))
if (flag_bc .gt. 0) then
allocate (C_bc(nlay))
allocate (N_bc(nlay))
C_bc = 0.
N_bc = 0.
endif
do i=1,nspecies
slit(i)%C_opm_frt = 0.
slit(i)%N_opm_frt = 0.
slit(i)%C_opm_crt = 0.
slit(i)%N_opm_crt = 0.
slit(i)%C_opm_tb = 0.
slit(i)%N_opm_tb = 0.
slit(i)%C_opm_stem = 0.
slit(i)%N_opm_stem = 0.
enddo
nlay2 = nlay+2
mfirst = 1
allocate (sh(mfirst:nlay2))
allocate (sv(mfirst:nlay2))
allocate (sb(mfirst:nlay2))
allocate (sbt(mfirst:nlay2))
allocate (t_cb(mfirst:nlay2))
allocate (t_cond(mfirst:nlay2))
allocate (h_cap(mfirst:nlay2))
allocate (sxx(mfirst:nlay2))
allocate (svv(mfirst:nlay2))
allocate (svva(mfirst:nlay2))
allocate (soh(mfirst:nlay2))
allocate (son(mfirst:nlay2+1))
help=1
C_opm = 0
allocate(fr_loss(nlay))
allocate(redis(nlay))
end subroutine alloc_soil
!------------------------------------------------------------------
subroutine stand_bal_ini
use data_stand
implicit none
integer i
allocate(diam_class(num_class, nspecies)); diam_class=0
allocate(diam_class_t(num_class, nspecies)); diam_class_t=0
allocate(diam_class_h(num_class,nspecies)); diam_class_h=0
allocate(diam_class_age(num_class,nspecies)); diam_class_age=0
allocate(diam_class_mvol(num_class,nspecies)); diam_class_mvol=0
allocate(diam_classm(num_class,nspecies)); diam_classm=0
allocate(diam_classm_h(num_class,nspecies)); diam_classm_h=0
allocate(height_class(num_class)); height_class =0
! array of potential litter (dead stems and twigs/branches for the next years
allocate(dead_wood(nspec_tree))
do i = 1,nspec_tree
allocate(dead_wood(i)%C_tb(lit_year))
allocate(dead_wood(i)%N_tb(lit_year))
allocate(dead_wood(i)%C_stem(lit_year))
allocate(dead_wood(i)%N_stem(lit_year))
dead_wood(i)%C_tb = 0.
dead_wood(i)%N_tb = 0.
dead_wood(i)%C_stem = 0.
dead_wood(i)%N_stem = 0.
enddo
end subroutine stand_bal_ini
!--------------------------------------------------------------
subroutine control_file ! saving simulation parameter and start conditions for each site
real buckdepth
character(8) actdate
character(10) acttime
character(150) site_help
integer help_ip, j
TYPE(Coh_Obj), Pointer :: help_coh ! pointer to cohort list
IF(site_nr==1) THEN
help_ip=site_nr
ELSE
help_ip=ip
END IF
! Write soil initialisation file
if (flag_mult8910) then
site_help = site_name1
else
site_help = site_name(help_ip)
endif
if (.not.flag_mult8910 .or. (flag_mult8910 .and. anh .eq. "1") .or. (flag_mult8910 .and. time_out .gt. 0.)) then
if (.not.flag_mult8910) then
WRITE (unit_soil,'(26A)') 'Layer',' Depth(cm)',' F-cap(mm)',' F-cap(Vol%)',' Wiltp(mm)', &
' Wiltp(Vol%)',' Pore vol.',' Skel.(Vol%)',' Density',' Spheat',' pH',' Wlam', &
' Water(mm)',' Water(Vol%)',' Soil-temp.',' C_opm g/m2', &
' C_hum g/m2',' N_opm g/m2',' N_hum g/m2',' NH4 g/m2',' NO3 g/m2',' humus part',' d_mass g/m2', ' Clay',' Silt',' Sand'
do i = 1,nlay
WRITE (unit_soil,'(I5,2F10.2,3F12.2,F10.2,F12.2,4F8.2,F10.2,F12.2, 5F11.2,2F9.4,2E12.4, 3F6.1)') i,depth(i),field_cap(i),f_cap_v(i),wilt_p(i), &
wilt_p_v(i),pv_v(i), skelv(i)*100., dens(i),spheat(i),phv(i),wlam(i), &
wats(i),watvol(i),temps(i),c_opm(i),c_hum(i),n_opm(i), n_hum(i),nh4(i),no3(i),humusv(i),dmass(i), clayv(i)*100., siltv(i)*100., sandv(i)*100.
end do
endif
! Write control file
call date_and_time(actdate, acttime)
unit_ctr = getunit()
open(unit_ctr,file=trim(dirout)//trim(site_help)//'.ctr'//anh,status='replace')
WRITE(unit_ctr,'(2A)') '*** Site name: ',site_name(help_ip)
WRITE(unit_ctr,'(2A)') ' Appendix ' ,anh
WRITE(unit_ctr,'(A,F7.2)') ' Longitude: ', long
WRITE(unit_ctr,'(A,F7.2)') ' Latitude: ', lat
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(10A)') ' ---- Version: v2.2 ---- '
WRITE(unit_ctr,'(10A)') ' Date: ',actdate(7:8),'.',actdate(5:6),'.',actdate(1:4), &
' Time: ',acttime(1:2),':',acttime(3:4)
WRITE(unit_ctr,'(A,A)') ' Simulation control file: ',trim(simfile)
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') '*** Data files:'
IF(flag_clim==1)then
WRITE(unit_ctr,'(A,A)') ' Climfile: ',trim(climfile(ip))
ELSE
WRITE(unit_ctr,'(A,A)') ' Climfile: ',trim(climfile(1))
endif
WRITE(unit_ctr,'(A,A)') ' Sitefile: ',trim(sitefile(help_ip))
WRITE(unit_ctr,'(A,A)') ' Start value file: ',trim(valfile(help_ip))
! Initialization of stand
IF( flag_multi==3 .OR. (site_nr>1 .AND. flag_stand>0) ) THEN
WRITE(unit_ctr,'(A,A)') ' Stand initialization: ',trim(treefile(ip))
ELSE IF( ip==1 .AND. flag_stand>0) THEN
WRITE(unit_ctr,'(A,A)') ' Stand initialization: ',trim(treefile(ip))
ELSE IF (flag_stand==0) THEN
WRITE(unit_ctr,'(A,A)') ' Stand initialization: none'
endif
IF (lmulti) WRITE(unit_ctr,'(A,A)') ' Stand identificator: ', adjustl(standid(ip))
WRITE(unit_ctr,*) ' '
IF(flag_mg.ne.0 .and. flag_mg.ne.5) then
WRITE(unit_ctr,'(A,A)') ' Management control file: ',trim(manfile(ip))
ELSE
WRITE(unit_ctr,'(A)') ' Management: none'
endif
WRITE(unit_ctr,'(A,A)') ' Deposition file: ',trim(depofile(ip))
WRITE(unit_ctr,'(A,A)') ' N reduction file: ',trim(redfile(ip))
WRITE(unit_ctr,'(A,A)') ' Litter initialisation file: ',trim(litfile(ip))
if (flag_stat .gt. 0) WRITE(unit_ctr,'(A,A)') ' File with measurements: ',trim(mesfile(1))
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') '*** Soil description '
WRITE(unit_ctr,'(A,I3)') ' Number of soil layers: ',nlay
WRITE(unit_ctr,'(A,I3)') ' Number of rooting layers: ',nroot_max
WRITE(unit_ctr,'(A,I3)') ' Ground water from layer: ',nlgrw
WRITE(unit_ctr,'(A,F5.1)') ' Evaporation depth (cm): ',w_ev_d
call bucket(bucks_100, bucks_root, buckdepth)
buckdepth = buckdepth/100
WRITE(unit_ctr,'(A,F5.2,A,F7.2)') ' Bucket size (mm), ', buckdepth,' m depth: ',bucks_100
WRITE(unit_ctr,'(A,F7.2)') ' Bucket size (mm) of rooting zone: ',bucks_root
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') '*** Soil water conditions'
WRITE(unit_ctr,'(12A)') 'Layer ','Depth(cm) ','F-cap(mm) ','F-cap(Vol%) ','Wiltp(mm) ', &
'Wiltp(Vol%) ','Pore vol. ','Density ','Spheat ','pH-value ',' Wlam',' skel. '
do i = 1,nlay
WRITE(unit_ctr,'(I5,12F10.2)') i,depth(i),field_cap(i),f_cap_v(i),wilt_p(i), &
wilt_p_v(i),pv_v(i),dens(i),spheat(i),phv(i),wlam(i),skelv(i)
end do
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') '*** Soil initial values'
WRITE(unit_ctr,'(9A)') 'Layer ','Water-cont. ','Soil-temp. ','C_opm ', &
'C_hum ','N_opm ','N_hum ','NH4-cont. ','NO3-cont '
do i=1,nlay
WRITE(unit_ctr,'(I5, 2F10.2, 6F10.4)') i,wats(i),temps(i),c_opm(i),c_hum(i),n_opm(i), &
n_hum(i),nh4(i),no3(i)
end do
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') ' N_tot C_tot N_antot N_humtot C_humtot C_opm_fol C_opm_tb C_opm_frt C_opm_crt C_opm_stem '
WRITE(unit_ctr,'(10F12.4)') N_tot, C_tot, N_an_tot, N_hum_tot, C_hum_tot, C_opm_fol, C_opm_tb, C_opm_frt, C_opm_crt, C_opm_stem
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)',advance='no') 'Mineralization constant of humus - humus layer (k_hum): '
WRITE(unit_ctr,'(F10.5)') k_hum
WRITE(unit_ctr,'(A)',advance='no') 'Mineralization constant of humus - mineral soil (k_hum_r): '
WRITE(unit_ctr,'(F10.5)') k_hum_r
WRITE(unit_ctr,'(A)',advance='no') 'Nitrification constant (k_nit): '
WRITE(unit_ctr,'(F10.5)') k_nit
WRITE(unit_ctr,*) ' '
if (flag_bc .gt.0) then
WRITE(unit_ctr,'(A)') '*** Biochar application '
WRITE(unit_ctr,'(A)') ' year C-content(%) C/N-ratio depth mass(kg/ha dry mass)'
do j = 1, n_appl_bc
WRITE(unit_ctr,'(I7,F14.1, F11.1, I7, F18.1)') &
y_bc(j), cpart_bc(j), cnv_bc(j), bc_appl_lay(j), C_bc_appl(j)
enddo
WRITE(unit_ctr,'(F10.5)')
endif
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') '*** Stand initialisation'
WRITE(unit_ctr,'(A)')' Coh x_fol x_frt x_sap x_hrt x_Ahb height x_hbole x_age n sp DC DBH'
help_coh => pt%first
DO WHILE (ASSOCIATED(help_coh))
WRITE(unit_ctr,'(I5,5f12.5,2f10.0,i7,f7.0,i7, 2f12.5)') help_coh%coh%ident, help_coh%coh%x_fol, help_coh%coh%x_frt, help_coh%coh%x_sap, help_coh%coh%x_hrt, &
help_coh%coh%x_Ahb, help_coh%coh%height, help_coh%coh%x_hbole, help_coh%coh%x_age, &
help_coh%coh%nTreeA,help_coh%coh%species, help_coh%coh%dcrb, help_coh%coh%diam
help_coh => help_coh%next
END DO
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') '*** Simulation control'
WRITE(unit_ctr,'(A66,I4)') 'Run option: ',flag_multi
WRITE(unit_ctr,'(A66,I4)') 'Start year: ',time_b
WRITE(unit_ctr,'(A66,I4)') 'Number of simulation years - year: ', year
WRITE(unit_ctr,'(A60,F12.1)') 'Patch size [m²] - kpatchsize: ',kpatchsize
WRITE(unit_ctr,'(A60,F12.1)') 'Thickness of leaf layers - dz: ',dz
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_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
WRITE(unit_ctr,'(A66,I4)') 'Disturbance (0-OFF, 1-ON ) - flag_dis: ',flag_dis
WRITE(unit_ctr,'(A66,I4)') 'Light absoption algorithm (1,2,3,4) - : ',flag_light
WRITE(unit_ctr,'(A66,I4)') 'Foliage-height relationship (0,1) - flag_folhei: ',flag_folhei
WRITE(unit_ctr,'(A66,I4)') 'Volume function trunc (0,1) - flag_volfunc: ',flag_volfunc
WRITE(unit_ctr,'(A66,I4)') 'Respiration model (0-0.5*NPP,1-organ specific) - flag_resp: ',flag_resp
WRITE(unit_ctr,'(A66,I4)') 'Limitation (0-NO,1-water, 2-N, 3-water+N) - flag_limi: ',flag_limi
WRITE(unit_ctr,'(A66,I4)') 'Flag for decomposition model - flag_decomp: ',flag_decomp
WRITE(unit_ctr,'(A66,I4)') 'Root spec. activity (0-const,1-varying) - flag_sign: ',flag_sign
WRITE(unit_ctr,'(A66,I4)') 'Water uptake function soil (1,2,3,4) - flag_wred: ',flag_wred
WRITE(unit_ctr,'(A66,I4)') 'Root distribution - flag_wurz: ',flag_wurz
WRITE(unit_ctr,'(A66,I4)') 'Heat conductance - flag_cond: ',flag_cond
WRITE(unit_ctr,'(A66,I4)') 'Interception - flag_int: ',flag_int
WRITE(unit_ctr,'(A66,I4)') 'Evapotranspiration - flag_eva: ',flag_eva
WRITE(unit_ctr,'(A66,I4)') 'CO2 (0-constant,1-historic increase,2-step change)- flag_co2: ',flag_co2
WRITE(unit_ctr,'(A66,I4)') 'Sort flag - flag_sort: ',flag_sort
WRITE(unit_ctr,'(A66,I4)') 'wpm flag - flag_wpm: ',flag_wpm
WRITE(unit_ctr,'(A66,I4)') 'Analysis of measurements - flag_stat: ',flag_stat
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A66,A)') 'Species parameter file: ',trim(specfile(help_ip))
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') '*** Species parameter description'
WRITE(unit_ctr,'(A51,I4)') ' Species number: ', nspecies
WRITE(unit_ctr,'(A51,I4)') ' Number of tree species: ', nspec_tree
WRITE(unit_ctr,*) ' ********** '
WRITE(unit_ctr,'(A25,A9,2X,A30)') 'Short Name', ' Spec-Nr', 'Latin Name '
WRITE(unit_ctr,*) ' '
do i=1,nspecies
WRITE(unit_ctr,'(A25,I9,2X,A30)') trim(spar(i)%species_short_name), i, spar(i)%species_name
enddo
WRITE(unit_ctr,*) ' ********** '
WRITE(unit_ctr,'(A51,15A16)') ' Species name: ', (trim(spar(i)%species_short_name),i=1,nspecies)
WRITE(unit_ctr,1010) ' Maximal age - max_age: ', (spar(i)%max_age,i=1,nspecies)
WRITE(unit_ctr,1010) ' Stress rec. time - yrec: ', (spar(i)%yrec,i=1,nspecies)
WRITE(unit_ctr,1010) ' Shade tolerance - stol: ', (spar(i)%stol,i=1,nspecies)
WRITE(unit_ctr,1000) ' Extinction coeff - pfext: ', (spar(i)%pfext,i=1,nspecies)
WRITE(unit_ctr,1000) ' Root activity rate - sigman: ', (spar(i)%sigman,i=1,nspecies)
WRITE(unit_ctr,1000) ' Respiration coeff - respcoeff: ', (spar(i)%respcoeff,i=1,nspecies)
WRITE(unit_ctr,1000) ' Growth resp. par. - prg: ', (spar(i)%prg,i=1,nspecies)
WRITE(unit_ctr,1000) ' Maint.resp.par./sapwood - prms: ', (spar(i)%prms,i=1,nspecies)
WRITE(unit_ctr,1000) ' Maint.resp.par./fineroot - prmr: ', (spar(i)%prmr,i=1,nspecies)
WRITE(unit_ctr,1000) ' Senesc.par. foliage - psf: ', (spar(i)%psf,i=1,nspecies)
WRITE(unit_ctr,1000) ' Senesc.par. sapwood - pss: ', (spar(i)%pss,i=1,nspecies)
WRITE(unit_ctr,1000) ' Senesc.par. fineroot - psr: ', (spar(i)%psr,i=1,nspecies)
WRITE(unit_ctr,1000) ' N/C ratio of biomass - pcnr: ', (spar(i)%pcnr,i=1,nspecies)
WRITE(unit_ctr,1000) ' N concentration of foliage - ncon_fol: ', (spar(i)%ncon_fol,i=1,nspecies)
WRITE(unit_ctr,1000) ' N concentration of fine roots - ncon_frt: ', (spar(i)%ncon_frt,i=1,nspecies)
WRITE(unit_ctr,1000) ' N concentration of coarse roots - ncon_crt: ', (spar(i)%ncon_crt,i=1,nspecies)
WRITE(unit_ctr,1000) ' N concentration of twigs and branches - ncon_tbc: ', (spar(i)%ncon_tbc,i=1,nspecies)
WRITE(unit_ctr,1000) ' N concentration of stemwood - ncon_stem: ', (spar(i)%ncon_stem,i=1,nspecies)
WRITE(unit_ctr,1000) ' Reallocation parameter of foliage - reallo_fol: ', (spar(i)%reallo_fol,i=1,nspecies)
WRITE(unit_ctr,1000) ' Reallocation parameter of fine root - reallo_frt: ', (spar(i)%reallo_frt,i=1,nspecies)
WRITE(unit_ctr,1000) ' Ratio of coarse wood - alphac: ', (spar(i)%alphac,i=1,nspecies)
WRITE(unit_ctr,1000) ' Coarse root fraction of coarse wood - cr_frac: ', (spar(i)%cr_frac,i=1,nspecies)
WRITE(unit_ctr,1000) ' Sapwood density - prhos: ', (spar(i)%prhos,i=1,nspecies)
WRITE(unit_ctr,1000) ' Proport.const.(pipe mod.) - pnus: ', (spar(i)%pnus,i=1,nspecies)
IF(flag_folhei==0) THEN
WRITE(unit_ctr,1000) ' Height growth parameter - pha: ', (spar(i)%pha,i=1,nspecies)
ELSEIF(flag_folhei==1) THEN
WRITE(unit_ctr,1000) ' Height growth par. 1 - pha_v1: ', (spar(i)%pha_v1,i=1,nspecies)
WRITE(unit_ctr,1000) ' Height growth par. 2 - pha_v2: ', (spar(i)%pha_v2,i=1,nspecies)
WRITE(unit_ctr,1000) ' Height growth par. 3 - pha_v3: ', (spar(i)%pha_v3,i=1,nspecies)
ELSE
WRITE(unit_ctr,'(A51,I3)') ' non valid flag value - flag_folhei : ',flag_folhei
ENDIF
WRITE(unit_ctr,1000) ' Height growth parameter coeff 1 - pha_coeff1: ', (spar(i)%pha_coeff1,i=1,nspecies)
WRITE(unit_ctr,1000) ' Height growth parameter coeff 2 - pha_coeff2: ', (spar(i)%pha_coeff2,i=1,nspecies)
WRITE(unit_ctr,1000) ' Crown radius - DBH ratio parameter a - crown_a: ', (spar(i)%crown_a,i=1,nspecies)
WRITE(unit_ctr,1000) ' Crown radius - DBH ratio parameter b - crown_b: ', (spar(i)%crown_b,i=1,nspecies)
WRITE(unit_ctr,1000) ' Crown radius - DBH ratio parameter c - crown_c: ', (spar(i)%crown_c,i=1,nspecies)
WRITE(unit_ctr,1000) ' Minimum specific leaf area - psla_min: ', (spar(i)%psla_min,i=1,nspecies)
WRITE(unit_ctr,1000) ' Light dep. specific leaf area - psla_a: ', (spar(i)%psla_a,i=1,nspecies)
WRITE(unit_ctr,1000) ' Efficiency parameter - phic: ', (spar(i)%phic,i=1,nspecies)
WRITE(unit_ctr,1000) ' N content - pnc: ', (spar(i)%pnc,i=1,nspecies)
WRITE(unit_ctr,1000) ' kco2_25: ', (spar(i)%kCO2_25,i=1,nspecies)
WRITE(unit_ctr,1000) ' ko2_25: ', (spar(i)%kO2_25,i=1,nspecies)
WRITE(unit_ctr,1000) ' CO2/O2 specif. value - pc_25: ', (spar(i)%pc_25,i=1,nspecies)
WRITE(unit_ctr,1000) ' Q10_kco2: ', (spar(i)%q10_kCO2,i=1,nspecies)
WRITE(unit_ctr,1000) ' Q10_ko2: ', (spar(i)%q10_kO2,i=1,nspecies)
WRITE(unit_ctr,1000) ' Q10_pc: ', (spar(i)%q10_pc,i=1,nspecies)
WRITE(unit_ctr,1000) ' Rd to Vm ratio - pb: ', (spar(i)%pb,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Inhibitor min temp. - PItmin: ', (spar(i)%PItmin,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Inhibitor opt temp. - PItopt: ', (spar(i)%PItopt,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Inhibitor max temp. - PItmax: ', (spar(i)%PItmax,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Inhibitor scaling factor - PIa: ', (spar(i)%PIa,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Promotor min temp. - PPtmin: ', (spar(i)%PPtmin,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Promotor opt temp. - PPtopt: ', (spar(i)%PPtopt,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Promotor max temp. - PPtmax: ', (spar(i)%PPtmax,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Promotor scaling factor - PPa: ', (spar(i)%PPa,i=1,nspecies)
WRITE(unit_ctr,1000) ' PIM: Promotor scaling factor - PPb: ', (spar(i)%PPb,i=1,nspecies)
WRITE(unit_ctr,1000) ' CSM: chilling base temp. - CSTbC: ', (spar(i)%CSTbC,i=1,nspecies)
WRITE(unit_ctr,1000) ' CSM: base temp. - CSTbT: ', (spar(i)%CSTbT,i=1,nspecies)
WRITE(unit_ctr,1000) ' CSM: scaling factor - CSa: ', (spar(i)%CSa,i=1,nspecies)
WRITE(unit_ctr,1000) ' CSM: scaling factor - CSb: ', (spar(i)%CSb,i=1,nspecies)
WRITE(unit_ctr,1000) ' TSM: base temp. - LTbT: ', (spar(i)%LTbT,i=1,nspecies)
WRITE(unit_ctr,1000) ' TSM: critical temperature sum - LTcrit: ', (spar(i)%LTcrit,i=1,nspecies)
WRITE(unit_ctr,1010) ' TSM: start day after 1.11. - Lstart: ', (spar(i)%Lstart,i=1,nspecies)
WRITE(unit_ctr,1000) ' usefd pheno model - Phmodel: ', (spar(i)%Phmodel,i=1,nspecies)
WRITE(unit_ctr,1000) ' End day for phenology - end_bb: ', (spar(i)%end_bb,i=1,nspecies)
WRITE(unit_ctr,1000) ' Fpar_mod - fpar_mod: ', (spar(i)%fpar_mod,i=1,nspecies)
WRITE(unit_ctr,1000) ' Intercep.cap. - ceppot_spec: ', (spar(i)%ceppot_spec,i=1,nspecies)
WRITE(unit_ctr,1000) ' photosynthesis response to N-limitation - Nresp: ', (spar(i)%Nresp,i=1,nspecies)
WRITE(unit_ctr,1000) ' Regeneration flag - regflag: ', (spar(i)%regflag,i=1,nspecies)
WRITE(unit_ctr,1000) ' Seedrate: ', (spar(i)%seedrate,i=1,nspecies)
WRITE(unit_ctr,1000) ' Seedmass: ', (spar(i)%seedmass,i=1,nspecies)
WRITE(unit_ctr,1000) ' Standard dev. of seedrate - seedsd: ', (spar(i)%seedsd,i=1,nspecies)
WRITE(unit_ctr,1000) ' all. parameter - seeda: ', (spar(i)%seeda,i=1,nspecies)
WRITE(unit_ctr,1000) ' all. parameter - seedb: ', (spar(i)%seedb,i=1,nspecies)
WRITE(unit_ctr,1000) ' all. parameter - pheight1: ', (spar(i)%pheight1,i=1,nspecies)
WRITE(unit_ctr,1000) ' all. parameter - pheight2: ', (spar(i)%pheight2,i=1,nspecies)
WRITE(unit_ctr,1000) ' all. parameter - pheight3: ', (spar(i)%pheight3,i=1,nspecies)
WRITE(unit_ctr,1000) ' all. parameter - pdiam1: ', (spar(i)%pdiam1,i=1,nspecies)
WRITE(unit_ctr,1000) ' all. parameter - pdiam2: ', (spar(i)%pdiam2,i=1,nspecies)
WRITE(unit_ctr,1000) ' all. parameter - pdiam3: ', (spar(i)%pdiam3,i=1,nspecies)
WRITE(unit_ctr,1000) ' decomp. parameter foliage - k_opm_fol: ', (spar(i)%k_opm_fol,i=1,nspecies)
WRITE(unit_ctr,1000) ' synth. parameter foliage - k_syn_fol: ', (spar(i)%k_syn_fol,i=1,nspecies)
WRITE(unit_ctr,1000) ' decomp. parameter fine roots - k_opm_frt: ', (spar(i)%k_opm_frt,i=1,nspecies)
WRITE(unit_ctr,1000) ' synth. parameter fine roots - k_syn_frt: ', (spar(i)%k_syn_frt,i=1,nspecies)
WRITE(unit_ctr,1000) ' decomp. parameter coarse roots - k_opm_crt: ', (spar(i)%k_opm_crt,i=1,nspecies)
WRITE(unit_ctr,1000) ' synth. parameter coarse roots - k_syn_crt: ', (spar(i)%k_syn_crt,i=1,nspecies)
WRITE(unit_ctr,1000) ' decomp. parameter twigs/branches - k_opm_tb: ', (spar(i)%k_opm_tb,i=1,nspecies)
WRITE(unit_ctr,1000) ' synth. parameter twigs/branches - k_syn_tb: ', (spar(i)%k_syn_tb,i=1,nspecies)
WRITE(unit_ctr,1000) ' decomp. parameter stem - k_opm_stem: ', (spar(i)%k_opm_stem,i=1,nspecies)
WRITE(unit_ctr,1000) ' synth. parameter dtem - k_syn_stem: ', (spar(i)%k_syn_stem,i=1,nspecies)
WRITE(unit_ctr,1000)
WRITE(unit_ctr,1000) ' spec_rl: ', (spar(i)%spec_rl,i=1,nspecies)
WRITE(unit_ctr,1000) ' tbase: ', (spar(i)%tbase,i=1,nspecies)
WRITE(unit_ctr,1000) ' topt: ', (spar(i)%topt,i=1,nspecies)
WRITE(unit_ctr,1000) ' bdmax_coef: ', (spar(i)%bdmax_coef,i=1,nspecies)
WRITE(unit_ctr,1000) ' porcrit_coef: ', (spar(i)%porcrit_coef,i=1,nspecies)
WRITE(unit_ctr,1000) ' ph_opt_max: ', (spar(i)%ph_opt_max,i=1,nspecies)
WRITE(unit_ctr,1000) ' ph_opt_min: ', (spar(i)%ph_opt_min,i=1,nspecies)
WRITE(unit_ctr,1000) ' ph_max: ', (spar(i)%ph_max,i=1,nspecies)
WRITE(unit_ctr,1000) ' ph_min : ', (spar(i)%ph_min ,i=1,nspecies)
WRITE(unit_ctr,1000) ' v_growth: ', (spar(i)%v_growth,i=1,nspecies)
WRITE(unit_ctr,1000)
WRITE(unit_ctr,1000) ' C/N ratio of foliage - cnr_fol: ', (spar(i)%cnr_fol,i=1,nspecies)
WRITE(unit_ctr,1000) ' C/N ratio of fine roots - cnr_frt: ', (spar(i)%cnr_frt,i=1,nspecies)
WRITE(unit_ctr,1000) ' C/N ratio of coarse roots - cnr_crt: ', (spar(i)%cnr_crt,i=1,nspecies)
WRITE(unit_ctr,1000) ' C/N ratio of twigs and branches - cnr_tbc: ', (spar(i)%cnr_tbc,i=1,nspecies)
WRITE(unit_ctr,1000) ' C/N ratio of stemwood - cnr_stem: ', (spar(i)%cnr_stem,i=1,nspecies)
WRITE(unit_ctr,1000)
WRITE(unit_ctr,1000) ' Reduction factor - RedN: ', (svar(i)%RedN, i=1,nspecies)
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,'(A)') '****** Model parameter ******'
WRITE(unit_ctr,1020) 'Optimum ratio of ci to ca [-] - Lambda: ',lambda
WRITE(unit_ctr,1020) 'Molar mass of carbon [g/mol] - Cmass: ',Cmass
WRITE(unit_ctr,1020) 'Minimum conductance [mol/(m2*d)] - gmin: ',gmin
WRITE(unit_ctr,1020) 'Shape of PS response curve - ps: ',ps
WRITE(unit_ctr,1020) 'Slope of N function at 20 °C [g(N) (mymol s-1)-1] - pn: ',pn
WRITE(unit_ctr,1020) 'Minimum N content [g/g] - nc0: ',nc0
WRITE(unit_ctr,1020) 'C3 quantum efficiency - qco2: ',qco2
WRITE(unit_ctr,1020) 'Scaling parameter - qco2a: ',qco2a
WRITE(unit_ctr,1020) 'Partial pressure of oxygen (kPa) - o2: ',o2
WRITE(unit_ctr,1020) 'Atmospheric CO2 content (mol/mol) - co2: ',co2_st
WRITE(unit_ctr,1020) 'Albedo of the canopy - pfref: ',pfref
WRITE(unit_ctr,1020) 'Part of C in biomass [-] - cpart: ',cpart
WRITE(unit_ctr,1020) 'Ratio of molecular weights of water and air - rmolw: ',rmolw
WRITE(unit_ctr,1020) 'Universal gas constant [J/mol/K] = [Pa/m3/K] - R_gas: ',R_gas
WRITE(unit_ctr,1020) 'von Karman''s constant [-] - c_karman: ',c_karman
WRITE(unit_ctr,1020) 'Specific heat of air at const. pressure [J/g/K] - c_air: ',c_air
WRITE(unit_ctr,1020) 'Psychrometer constant [hPa/K] - psycro: ',psycro
WRITE(unit_ctr,1020) 'Breast height for inventory measurements [cm] - h_breast: ',h_breast
WRITE(unit_ctr,1020) 'Height for sapling allometry - h_sapini: ',h_sapini
WRITE(unit_ctr,1020) 'Min. diff. b. height of crown base and breast height- h_bo_br_diff: ',h_bo_br_diff
WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p1_co2: ',p1_co2
WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p2_co2: ',p2_co2
WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p3_co2: ',p3_co2
WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p4_co2: ',p4_co2
WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p5_co2: ',p5_co2
WRITE(unit_ctr,1020) 'Parameter variable for calculation of historical CO2 scenario - p1_co2h: ',p1_co2h
WRITE(unit_ctr,1020) 'Parameter variable for calculation of historical CO2 scenario - p2_co2h: ',p2_co2h
WRITE(unit_ctr,1020) 'Parameter variable for calculation of historical CO2 scenario - p3_co2h: ',p3_co2h
WRITE(unit_ctr,1020) 'Parameter variable for calculation of historical CO2 scenario - p4_co2h: ',p4_co2h
WRITE(unit_ctr,1020) 'Threshold of air temperature for snow accumulation [°C] - temp_snow: ',temp_snow
WRITE(unit_ctr,1020) 'Parameter for calculation of transpiration demand - alfm: ',alfm
WRITE(unit_ctr,1020) 'Parameter for calculation of transpiration demand [mol/(m2*d)] - gpmax: ',gpmax
WRITE(unit_ctr,1020) 'Parameter for growing degree day calculation - thr_gdd: ',thr_gdd
IF (flag_multi==2) THEN
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,*) 'runs with climate scenarios produced by adding summands to every daily temperature'
WRITE(unit_ctr,*) 'and modifying every single precipitation value by a multiplicand'
WRITE(unit_ctr,*) 'run ident deltaT delta P factor'
ENDIF
! mangament parameter adaptation management
IF (flag_mg.eq.2. .and. flag_reg .eq. 0) then
WRITE(unit_ctr,*) ' '
WRITE(unit_ctr,*) '***Managment parameter case flag_mg = 2 (user specified) ***'
WRITE(unit_ctr,'(A35,4F15.5)') 'height for management control(cm)', ho1,ho2,ho3,ho4
WRITE(unit_ctr,'(A35,5I15)') 'management flags thr1-thr5' , thr1,thr2, thr3,thr4,thr5
WRITE(unit_ctr,'(A35,F15.5)') 'height for directional felling', thr6
WRITE(unit_ctr,'(A35,I15)') 'measure at rotation', thr7
WRITE(unit_ctr,'(A35,I15)') 'regeneration measure', mgreg
WRITE(unit_ctr,'(A35,F15.5)') 'lower/upper limit of height(cm)', limit
WRITE(unit_ctr,'(A35,I15)') 'number of years between thinning',thinstep
WRITE(unit_ctr,'(A35,F15.5)') 'rel. value for directional felling', direcfel
WRITE(unit_ctr,'(A35,5F15.5)')'number of Zielbaeume(spec.)', (zbnr(i),i=1,nspec_tree)
WRITE(unit_ctr,'(A35,5F15.5)')'rel. value for tending of pl.',(tend(i), i =1,nspec_tree)
WRITE(unit_ctr,'(A35,5I15)')'rotation ',(rot(i), i =1,nspec_tree)
WRITE(unit_ctr,'(A35,5I15)')'age of nat./pl. regeneration',(regage(i), i =1,nspec_tree)
end IF
IF (flag_multi .ne. 2.and. flag_mg.ne.2 .and. flag_reg .eq.0) close(unit_ctr)
endif ! flag_mult8910
1000 FORMAT (A51,15 F16.5)
1010 FORMAT (A51,15 I16)
1020 FORMAT(A70,F15.5)
end subroutine control_file
end subroutine prepare_site
!******************************************************************************
SUBROUTINE readbudb
use data_simul
use data_species
use data_stand
implicit none
DO ns=1,nspecies
IF(spar(ns)%phmodel==4) THEN
WRITE(*,*) 'Please type the day of budburst for 4C species number ',ns,':'
READ(*,*) svar(ns)%ext_daybb
ENDIF
ENDDO
END subroutine readbudb
!******************************************************************************
SUBROUTINE readdepo
use data_climate
use data_depo
use data_out
use data_simul
use data_site
implicit none
character text
integer hx, unit_dep, i,j,ios, ii
!integer realrec
integer id,im,iy,itz1, itz2, hyear1, hyear2, hyear3, hy
logical ex
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
NOd = NOdep(ip) ! concentration mg/l
NHd = NHdep(ip) ! concentration mg/l
return
endif
NOd = 0.
NHd = 0.
if (.not.flag_mult8910) print *
inquire (File = depofile(ip), exist = ex) ! test whether file exists
IF(ex .eqv. .false.) then
if (.not.flag_mult8910) then
hx = 0
print *,' >>>FORESEE message: Cannot find deposition data - all data set to zero!'
CALL error_mess(hx,'Cannot find deposition data - all data set to ',REAL(hx))
endif
else
if (.not.flag_mult8910) print *, ' >>>FORESEE message: Now reading DEPOSITION data from file, please wait...'
! now read data from file
unit_dep = getunit()
OPEN (unit_dep,FILE=depofile(ip),IOSTAT=ios,STATUS='OLD',ACTION='READ')
flag_depo = 1
read(unit_dep,*) text
select case (text)
case ('C', 'c') ! concentrations mg/l
flag_depo = 2
read(unit_dep,*) text
case ('Y', 'y') ! Yearly constant deposition mg/m2
flag_depo = 3
read(unit_dep,*) text
case ('A', 'a') ! Annual sum of deposition g/m2
flag_depo = 4
read(unit_dep,*) text
end select
do
IF (text .ne. '!') then
backspace(unit_dep)
exit
endif
read(unit_dep,*) text
enddo
! assignment of dates
! fill in missing values with current values until current date
! fill in missing values at the end
hyear1 = 0
hyear2 = 0
hyear3 = 1
itz1 = 1
itz2 = 1
select case (flag_depo)
case(4)
do while ((ios .eq. 0) .and. (hyear1 .lt. year))
read(unit_dep,*,iostat=ios) iy, hNO, hNH
if (ios .eq.0) then
if (iy .gt. time_b+year) then
hyear1 = year
else
hyear1 = iy - time_b + 1
endif
if ((hyear1 .le. year) .and. (hyear1 .gt. 0)) then ! save from simulation start year onwards
do i = 1,366
NOd(i,hyear1) = hNO * 1000./366. ! report of year [g/m2] in daily values [mg/m2]
NHd(i,hyear1) = hNH * 1000./366.
enddo
hy = hyear1-1
do while ((hy .gt. hyear2) .and. (hy .gt. 0))
do i = 366, 1, -1
NOd(i,hy) = hNO * 1000./366.
NHd(i,hy) = hNH * 1000./366.
enddo
hy = hy - 1
enddo
hyear2 = hyear1
endif ! 0 < hyear1 < year
else ! ios .ne. 0
if (hyear1 .le. 0) then
hyear1 = 1
hyear2 = 1
endif
continue
endif ! ios = 0
enddo
case default
do while ((ios .eq. 0) .and. (hyear1 .lt. year))
read(unit_dep,*,iostat=ios) id,im,iy, hNO, hNH
if (ios .eq.0) then
call daintz(id,im,iy,itz1)
if (iy .gt. time_b+year) then
hyear1 = year
else
hyear1 = iy - time_b + 1
endif
if ((hyear1 .le. year) .and. (hyear1 .gt. 0)) then ! save from simulation start year onwards
NOd(itz1,hyear1) = hNO
NHd(itz1,hyear1) = hNH
select case (flag_depo)
case(1,2)
if (hyear1 .eq. hyear3) then
if (itz1 .gt. 1) then
do i = itz1-1, itz2, -1
NOd(i,hyear1) = hNO
NHd(i,hyear1) = hNH
enddo
endif
else
if (itz2 .lt. recs(hyear3)) then
do i = itz2+1, recs(hyear3)
NOd(i,hyear3) = hNO
NHd(i,hyear3) = hNH
enddo
endif
itz2 = 1
if (itz1 .gt. 1) then
do i = itz1-1, itz2, -1
NOd(i,hyear1) = hNO
NHd(i,hyear1) = hNH
enddo
endif
hy = hyear1-1
do while ((hy .gt. hyear3) .and. (hy .gt. 0))
do i = 366, 1, -1
NOd(i,hy) = hNO
NHd(i,hy) = hNH
enddo
hy = hy - 1
enddo
endif ! hyear1 .eq. hyear3
hyear3 = hyear1
itz2 = itz1
hyear2 = hyear3
case(3) ! fill in of constant year values
do i = 1,366
NOd(i,hyear1) = hNO
NHd(i,hyear1) = hNH
enddo
hy = hyear1-1
do while ((hy .gt. hyear2) .and. (hy .gt. 0))
do i = 366, 1, -1
NOd(i,hy) = hNO
NHd(i,hy) = hNH
enddo
hy = hy - 1
enddo
hyear2 = hyear1
itz2 = 366
end select ! flag_depo 1-3
endif ! 0 < hyear1 < year
else ! ios .ne. 0
if (hyear1 .le. 0) then
hyear1 = 1
hyear2 = 1
endif
continue
endif ! ios = 0
enddo
end select ! flag_depo
! fill in of the missing data at the end
select case (flag_depo)
case (3)
if (hyear1 .lt. year) then
hy = hyear1+1
do while (hy .le. year)
do i = 366, 1, -1
NOd(i,hy) = hNO
NHd(i,hy) = hNH
enddo
hy = hy + 1
enddo
else ! if date is outside the simulation period, it will be completly filled in
do j = 1, year
do i = 1, 366
NOd(i,j) = hNO
NHd(i,j) = hNH
enddo
enddo
endif
case default
if (hyear2 .le. year) then
if (itz2 .lt. recs(hyear2)) then
if (.not.flag_mult8910) then
hx = iy
CALL error_mess(hx,' Not enough data records in deposition file, iostat = ',REAL(ios))
WRITE (unit_err,*) ' >>>FORESEE message: Fill next values with same data '
WRITE (unit_err,'(A,2I4,A,2I4)')' from internal simulation time', itz2, hyear2, ' to', recs(hyear2), year
endif
do j = hyear2, year
ii = 1
if (j .eq. hyear2) ii = itz2
do i = ii, 366
NOd(i,j) = hNO
NHd(i,j) = hNH
enddo
enddo
else
hy = hyear2+1
do while (hy .le. year)
do i = 366, 1, -1
NOd(i,hy) = hNO
NHd(i,hy) = hNH
enddo
hy = hy + 1
enddo
endif
else ! if date is outside the simulation period, it will be completly filled in
do j = 1, year
do i = 1, 366
NOd(i,j) = hNO
NHd(i,j) = hNH
enddo
enddo
endif
end select
close (unit_dep)
endif
write (*,*)
END subroutine readdepo
!******************************************************************************
SUBROUTINE readredN
use data_out
use data_site
use data_species
use data_stand
use data_simul
implicit none
character text
integer hx, unit_red, i,ios
logical ex
if (.not.flag_mult8910) print *
if (flag_multi .lt. 8) then
inquire (File = "./input/.", exist = ex) ! test whether file exists
inquire (File = redfile(ip), exist = ex) ! test whether file exists
IF(ex .eqv. .false.) then
print *,' >>>FORESEE message: Cannot find data of RedN - internal calculation'
hx = 0
CALL error_mess(hx,'Cannot find data of RedN - internal calculation ',REAL(hx))
else
print *, ' >>>FORESEE message: Now reading RedN data from file, please wait...'
unit_red = getunit()
OPEN (unit_red,FILE=redfile(ip),IOSTAT=ios,STATUS='OLD',ACTION='READ')
DO
READ(unit_red,*) text
IF (text .ne. '!') THEN
backspace(unit_red)
EXIT
ENDIF
ENDDO
read (unit_red,*,iostat=ios) (svar(i)%RedN, i=1,nspecies)
close (unit_red)
endif ! ex
else
do i = 1, nspecies
svar(i)%RedN = RedN_list(i, ip)
enddo
endif ! flag_multi
IF(flag_limi==0 .OR. flag_limi==1) THEN
DO i=1,nspecies
svar(i)%RedN = 1.
ENDDO
ENDIF
do i = 1,nspecies
if (svar(i)%RedN .lt. 0) then ! no values; internal calculation
if (flag_multi .lt. 8) then
print *,' >>>FORESEE message: Cannot find data of RedN - internal calculation for', spar(i)%species_short_name
write (unit_err, '(A,I3,1X,A)') 'Cannot find data of RedN - internal calculation for species ',i, spar(i)%species_short_name
endif
flag_redn = .TRUE.
endif
enddo
if (.not.flag_mult8910) write (*,*)
END subroutine readredN
!******************************************************************************
SUBROUTINE readlit
!use data_climate
use data_out
use data_soil_cn
use data_species
use data_stand
use data_simul
implicit none
character text
integer unit_lit, i,ios
integer nspec_lit
logical ex
real help, hx
real, dimension(22) :: helpin
flag_lit = 0
if (flag_mult8910) then
inquire (File = litfile(1), exist = ex) ! test whether file exists
else
print *
inquire (File = litfile(ip), exist = ex) ! test whether file exists
endif
IF(ex .eqv. .false.) then
if (.not.flag_mult8910) then
print *,' >>>FORESEE message: Cannot find data of litter initialisation - internal calculation'
hx = 0.
write (unit_err,*)
write (unit_err,*) 'Cannot find data of litter initialisation - internal calculation '
endif
else
if (.not.flag_mult8910) print *, ' >>>FORESEE message: Now reading litter initialisation data from file, please wait...'
! now read data from file
unit_lit = getunit()
OPEN (unit_lit,FILE=litfile(ip),IOSTAT=ios,STATUS='OLD',ACTION='READ')
do
read(unit_lit,*) text
IF (text .ne. '!') then
backspace(unit_lit)
exit
endif
enddo
helpin = 0.
slit%C_opm_fol = 0.
read (unit_lit,*) nspec_lit
read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_fol, i=1,nspec_lit)
read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_tb , i=1,nspec_lit)
read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_frt(1), i=1,nspec_lit)
read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_crt(1), i=1,nspec_lit)
read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_stem,i=1,nspec_lit)
flag_lit = 1
help = 0.
hx = 1.
do i=1,nspecies
if (slit(i)%C_opm_fol .gt. 0) then
totfol_lit = totfol_lit + slit(i)%C_opm_fol
totfrt_lit = totfrt_lit + slit(i)%C_opm_frt(1)
tottb_lit = tottb_lit + slit(i)%C_opm_tb
totcrt_lit = totcrt_lit + slit(i)%C_opm_crt(1)
totstem_lit = totstem_lit + slit(i)%C_opm_stem
else
hx = -1.
endif
enddo
help = totfol_lit
if ((help .gt. 0.) .or. (hx .gt. 0) .and. .not.flag_mult8910) then
CALL error_mess(0,'Using data of litter initialisation from file '//trim(litfile(ip)), hx)
else
! no values; internal calculation of litter initialisation
if (.not.flag_mult8910) then
print *,' >>>FORESEE message: No data of litter initialisation - internal calculation'
hx = 0.
CALL error_mess(0,'No data of litter initialisation - internal calculation ', hx)
endif
flag_lit = 0
endif
close (unit_lit)
endif ! ex
if (.not.flag_mult8910) write (*,*)
END subroutine readlit
!******************************************************************************
subroutine prepare_climate
! read climate file
use data_climate
use data_out
use data_simul
use data_stand
implicit none
type clifile ! new data type for all climate parameters
integer :: day,mon,ye
real :: m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11
end type clifile
type (clifile), allocatable,dimension(:,:) :: climall !variable for data type climfile
character(1) c
character :: text
integer :: i,j,ios, unit_cli
integer :: realrec = 0
integer :: repflag = 0
logical :: ex
if (.not.flag_mult8910) then
print *, ' '
print *, ' Input of climate data: '
endif
call testfile(climfile(ip),ex) !input filename, test whether file exists
IF(ex .eqv. .false.) then
print *,' >>>FORESEE message: Cannot find climate data - program STOP!'
stop
endif
if (.not.flag_mult8910) print *, ' >>>FORESEE message: Now reading CLIMATE data from file, please wait...'
!now read data from file
unit_cli = getunit()
OPEN (unit_cli,FILE=climfile(ip),IOSTAT=ios,STATUS='OLD',ACTION='READ')
allocate (recs (1:year))
allocate (dd (1:366,1:year));allocate (mm (1:366, 1:year))
allocate (yy (1:year));allocate (tp (-2:366,1:year))
allocate (hm (0:366,1:year));allocate (prc (0:366,1:year))
allocate (prs (0:366,1:year));allocate (rd (0:366,1:year))
allocate (tn (0:366,1:year))
allocate (tx (0:366,1:year))
allocate (vp (0:366,1:year))
allocate (sdu (0:366,1:year))
allocate (wd (0:366,1:year))
allocate (sde (0:366,1:year))
allocate (bw (0:366,1:year))
dd = -99.9
mm = -99.9
yy = -99.9
tn = -99.9
tx = -99.9
wd = -99.9 ! wind initialisation
IF (index(climfile(ip),'.cli') .ne. 0) then
flag_climtyp = 1
do
read(unit_cli,*) text
IF (text .ne. '!') then
IF (text .eq. 'N') then
flag_climtyp = 2
else IF(text.eq.'T') then
flag_climtyp = 3
else
backspace(unit_cli)
exit
endif
endif
enddo
else if (index(climfile(ip),'.txt') .ne. 0) then
flag_climtyp = 4
else
flag_climtyp = 5
end IF
call read_cli
close(unit_cli)
if (flag_end .gt. 0) return
IF (realrec < year .and. repflag == 0) then
year = realrec
else
IF (repflag == 1) then
call climfill
end IF
end IF
med_rad1 = 0.
do j = 1, year-1
tp(0,j+1) = tp(recs(j),j)
tp(-1,j+1)= tp(recs(j)-1,j)
tp(-2,j+1)= tp(recs(j)-2,j)
hm(0,j+1) = hm(recs(j),j);prc(0,j+1) = prc(recs(j),j);prs(0,j+1) = prs(recs(j),j)
rd(0,j+1) = rd(recs(j),j)
wd(0,j+1) = wd(recs(j),j)
bw(0,j+1) = bw(recs(j),j)
vp(0,j+1) = vp(recs(j),j)
sdu(0,j+1) = sdu(recs(j),j)
sde(0,j+1) = sde(recs(j),j)
tx(0,j+1) = tx(recs(j),j)
tn(0,j+1) = tn(recs(j),j)
if( yy(j) .eq.time_b) then
do i=1, recs(j)
med_rad1 = med_rad1 + rd(i, j)
end do
med_rad1 = med_rad1/recs(1)
end if
end do
tp(-2,1) = tp(1,1); tp(-1,1) = tp(1,1); tp(0,1) = tp(1,1)
hm(0,1) = hm(1,1);prc(0,1) = prc(1,1);prs(0,1) = prs(1,1)
rd(0,1) = rd(1,1)
wd(0,1)=wd(1,1)
vp(0,1) = vp(1,1)
bw(0,1) = bw(1,1)
tn(0,1) = tn(1,1)
tx(0,1) = tx(1,1)
sdu(0,1) =sdu(1,1)
sde(0,1) = sde(1,1)
contains
!--------------------------------------------------------------
subroutine read_dwd
character(3) text
integer help, help1, help2, help3
allocate (climall (0:366,1:year))
j=1
c = 'n'
do
IF (j > year) then
realrec = year
exit
end IF
if (.not.flag_mult8910) print *, 'Year ',j
read(unit_cli,*) text
if(text.ne.'ta ') then
backspace(unit_cli)
end if
do i = 1, 366
read (unit_cli,*,IOSTAT=ios) climall(i,j)
help2 = climall(i,j)%day
help3 = climall(i,j)%mon
help = climall(i,j)%ye
help1 = climall(i-1,j)%ye
if (help.eq.2099 .and.help1.eq.2100.and. i.eq.366) then
end if
end do
IF (climall(365,j)%ye == climall(366,j)%ye) then
recs(j) = 366
else
backspace unit_cli
climall(366,j)%day = 0
climall(366,j)%mon = 0
climall(366,j)%ye = 0
recs(j) = 365
help = help-1
end IF
IF (j < year .and. ios < 0 .and. c .eq. 'n') then
realrec = j
if (.not.flag_mult8910) then
print *, ' >>>FORESEE message: Not enough climate data records in file!'
call error_mess(0,'read_cli: Not enough data records in climate file; number of complete years: ',real(realrec))
write(unit_err,'(A,I5)')' read_cli: Fill next values with same from first year, day: ',i_exit
write(unit_err,'(A,I5)')' read_cli: Fill next values with same data up to years: ',year
repflag = 1
exit
endif
else if(j.eq.year.and.ios < 0) then
realrec = year
exit
end IF
j=j+1
if(help.lt.time_b) j = j-1
end do
do j = 1, realrec
yy(j) = climall(1,j)%ye
do i = 1, recs(j)
dd(i,j) = climall(i,j)%day
mm(i,j) = climall(i,j)%mon
tx(i,j) = climall(i,j)%m1
tp(i,j) = climall(i,j)%m2
tn(i,j) = climall(i,j)%m3
prc(i,j) = climall(i,j)%m4
hm(i,j) = climall(i,j)%m5
prs(i,j) = climall(i,j)%m6
vp(i,j) = climall(i,j)%m7
sdu(i,j) = climall(i,j)%m8
bw(i,j) = climall(i,j)%m9
rd(i,j) = climall(i,j)%m10
wd(i,j) = climall(i,j)%m11
end do
end do
close(9)
deallocate (climall)
end subroutine read_dwd
!--------------------------------------------------------------
subroutine read_cli
implicit none
integer :: testtext, hp
character(11) :: text2
character(4) :: text
testtext=0
c = 'n'
j = 1
hp = 0
read(unit_cli,'(A)') text2
hp = index(text2,'.')
backspace(unit_cli)
do
IF(j > year) exit
select case(flag_climtyp)
case (1)
do i=1,366
if (hp .gt. 0) then
read(unit_cli,*,iostat=ios) text2,tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j)
text = text2(1:2)
write (text,'(A)') text2(1:2)
read (text,*) dd(i,j)
write (text,'(A)') text2(4:5)
read (text,*) mm(i,j)
write (text,'(A)') text2(7:10)
read (text,*) yy(j)
else
read(unit_cli,*,iostat=ios) dd(i,j),mm(i,j),yy(j),tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j)
endif ! hp
i_exit = i
if ((dd(i,j) .eq. 31) .and. (mm(i,j) .eq. 12)) then
recs(j) = i
write (*,*) 'Year ',j, yy(j)
realrec = j
if (j .eq. year) ios = -10
exit
endif
if (ios .ne. 0) exit
end do
case (2)
do i=1,366
read(unit_cli,*) dd(i,j),mm(i,j),yy(j),&
tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j),wd(i,j)
i_exit = i
if ((dd(i,j) .eq. 31) .and. (mm(i,j) .eq. 12)) then
recs(j) = i
write (*,*) 'Year ',j, yy(j)
realrec = j
if (j .eq. year) ios = -10
exit
endif
if (ios .ne. 0) exit
end do
case (3)
do i=1,366
if (hp .gt. 0) then
read(unit_cli,*,iostat=ios) text2, &
tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j),wd(i,j), tx(i,j),tn(i,j)
text = text2(1:2)
write (text,'(A)') text2(1:2)
read (text,*) dd(i,j)
write (text,'(A)') text2(4:5)
read (text,*) mm(i,j)
write (text,'(A)') text2(7:10)
read (text,*) yy(j)
else
read(unit_cli,*,iostat=ios) dd(i,j),mm(i,j),yy(j),&
tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j),wd(i,j), tx(i,j),tn(i,j)
endif
i_exit = i
if ((dd(i,j) .eq. 31) .and. (mm(i,j) .eq. 12)) then
recs(j) = i
write (*,*) 'Year ',j, yy(j)
realrec = j
if (j .eq. year) ios = -10
exit
endif
if (ios .ne. 0) exit
end do
case (4) ! suffix 'txt'
if (j .eq. 1 .and. testtext.eq.0) then
read(unit_cli,*) text
testtext = 1
end if
do i=1,366
read(unit_cli,*,iostat=ios) dd(i,j),mm(i,j),yy(j),&
tx(i,j),tp(i,j),tn(i,j),prc(i,j),hm(i,j),prs(i,j),rd(i,j),wd(i,j)
i_exit = i
if ((dd(i,j) .eq. 31) .and. (mm(i,j) .eq. 12)) then
recs(j) = i
write (*,*) 'Year ',j, yy(j)
realrec = j
if (j .eq. year) ios = -10
exit
endif
if (ios .ne. 0) exit
end do
case (5 )
call read_dwd
exit
end select
IF (realrec .lt. year .and. ios .ne. 0 .and. c .eq. 'n') then
if (dd(i_exit,j) .gt. 0) i_exit = i_exit+1
if (i_exit .ge. 365) i_exit = 1
repflag = 1
if (.not.flag_mult8910) then
print *, ' >>>FORESEE message: Not enough data records in file'
print *, ' IOSTAT = ', ios
WRITE (*,'(A,I5)') ' >>>FORESEE message: Fill next values with same data from day number', i_exit
CALL error_mess(0,'read_cli: Not enough data records in meteorology file; number of complete years: ',real(realrec))
write(unit_err,'(A,I5)')' read_cli: Fill next values with same from first year, day: ',i_exit
write(unit_err,'(A,I5)')' read_cli: Fill next values with same data up to years: ',year
exit
endif
end if
if (ios .ne. 0) exit
if (yy(j) .ge. time_b) then
if ((j .eq. 1) .and. (yy(j) .gt. time_b)) then
CALL error_mess(0,'read_cli: No climate data in meteorology file for year ',real(time_b))
flag_end = 6
return
endif
j = j+1
endif
end do
end subroutine read_cli
!--------------------------------------------------------------
subroutine climfill
integer istart
istart = i_exit
if(istart.eq.0) istart =istart +1
do j=realrec+1,year
print *,'Year ',j
yy(j)=yy(j-realrec)
recs(j)=recs(j-realrec)
do i=istart,366
dd(i,j) = dd(i,j-realrec)
mm(i,j) = mm(i,j-realrec)
tp(i,j) = tp(i,j-realrec)
hm(i,j) = hm(i,j-realrec)
prc(i,j) = prc(i,j-realrec)
prs(i,j) = prs(i,j-realrec)
rd(i,j) = rd(i,j-realrec)
wd(i,j) = wd(i,j-realrec)
tx(i,j) = tx(i,j-realrec)
tn(i,j) = tn(i,j-realrec)
end do
end do
end subroutine climfill
END subroutine prepare_climate
!**************************************************************
SUBROUTINE store_para(hpara, simpara, parerr)
use data_simul
use data_out
use data_par
use data_species
use data_soil_cn
use data_stand
use data_tsort
implicit none
integer inum
real hpara, parerr
character(100):: simpara, hchar1
integer, external :: array_num
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' store_para'
parerr = 0.
if (trim(simpara) .eq. 'year') then
year=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'time_b') then
time_b=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'kpatchsize') then
kpatchsize=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'dz') then
dz=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'ns_pro') then
ns_pro=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_mort') then
flag_mort=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_reg') then
flag_reg=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_stand') then
flag_stand=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_sveg') then
flag_sveg=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_mg') then
flag_mg=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_dis') then
flag_dis=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_light') then
flag_light=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_folhei') then
flag_folhei=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_volfunc') then
flag_volfunc=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_resp') then
flag_resp=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_limi') then
flag_limi=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_sign') then
flag_sign=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_decomp') then
flag_decomp=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_wred') then
flag_wred=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_wurz') then
flag_wurz=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_cond') then
flag_cond=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_int') then
flag_int=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_eva') then
flag_eva=hpara
parerr = 1.
return
endif
if ((trim(simpara) .eq. 'flag_co2') .or.(trim(simpara) .eq. 'flag_CO2')) then
flag_co2=hpara
parerr = 1.
return
endif
if (adjustl(trim(simpara)) .eq. 'flag_sort') then
flag_sort = hpara
parerr = 1.
return
endif
if (adjustl(trim(simpara)) .eq. 'flag_wpm') then
flag_wpm = hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'time_out') then
time_out=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_dayout') then
flag_dayout=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_cohout') then
flag_cohout=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'flag_sum') then
flag_sum=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'k_hum') then
k_hum=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'k_hum_r') then
k_hum_r=hpara
parerr = 1.
return
endif
if (trim(simpara) .eq. 'k_nit') then
k_nit=hpara
parerr = 1.
return
endif
if (adjustl(trim(simpara)) .eq. 'alfm') then
alfm = hpara
parerr = 1.
return
endif
if (adjustl(trim(simpara)) .eq. 'gpmax') then
gpmax = hpara
parerr = 1.
return
endif
if (adjustl(trim(simpara)) .eq. 'alfm') then
alfm = hpara
parerr = 1.
return
endif
! Species parameter
hchar1 = adjustl(simpara)
inum = array_num(hchar1)
if (hchar1(1:9) .eq. 'k_opm_fol') then
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%k_opm_fol = hpara
parerr = 1.
return
endif
endif
if (hchar1(1:9) .eq. 'k_opm_frt') then
inum = array_num(hchar1)
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%k_opm_frt = hpara
parerr = 1.
return
endif
endif
if (hchar1(1:9) .eq. 'k_syn_fol') then
inum = array_num(hchar1)
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%k_syn_fol = hpara
parerr = 1.
return
endif
endif
if (hchar1(1:9) .eq. 'k_syn_frt') then
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%k_syn_frt = hpara
parerr = 1.
return
endif
endif
if (hchar1(1:3) .eq. 'psf') then
inum = array_num(hchar1)
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%psf = hpara
parerr = 1.
return
endif
endif
if (hchar1(1:7) .eq. 'Phmodel') then
inum = array_num(hchar1)
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%Phmodel = hpara
parerr = 1.
return
endif
endif
if ((hchar1(1:4) .eq. 'pnus') .or. (hchar1(1:4) .eq. 'Pnus')) then
inum = array_num(hchar1)
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%pnus = hpara
parerr = 1.
return
endif
endif
if ((hchar1(1:4) .eq. 'RedN') .or. (hchar1(1:4) .eq. 'redn')) then
inum = array_num(hchar1)
if (inum .gt. 0 .and. inum .le. nspecies) then
svar(inum)%RedN = hpara
parerr = 1.
return
endif
endif
if (hchar1(1:4) .eq. 'prms') then
inum = array_num(hchar1)
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%prms = hpara
parerr = 1.
return
endif
endif
if (hchar1(1:4) .eq. 'prmr') then
inum = array_num(hchar1)
if (inum .gt. 0 .and. inum .le. nspecies) then
spar(inum)%prmr = hpara
parerr = 1.
return
endif
endif
END subroutine store_para
!**************************************************************
integer FUNCTION array_num(string)
! reads the field numbre out of an array and hands it back as integer
implicit none
integer ipos1, ipos2, inum
character (100) string
character (10) help, hchar
ipos1 = scan(string, '(' )
ipos2 = scan(string, ')' )
ipos1 = ipos1+1
ipos2 = ipos2-1
hchar = string(ipos1:ipos2)
write(help,'(A3)') hchar
read(help,*) inum
array_num = inum
end function array_num