Forked from
4C / FORESEE
182 commits behind the upstream repository.
-
Petra Lasch-Born authoredPetra Lasch-Born authored
simul.f 13.11 KiB
!*****************************************************************!
!* *!
!* 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
logical lhelp
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
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)? '
read *,a
IF (a .eq. 'n' .or. a .eq. 'N') then
flag_dayout = 0
ENDIF
endif ! lhelpENDIF
! 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 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) 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) CALL dis_manag
! 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).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
!**************************************************************