source_code/version2.2_windows/amod_soil.f

@@ -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
 

source_code/version2.2_windows/day_ini.f

@@ -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

source_code/version2.2_windows/initia.f

@@ -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

source_code/version2.2_windows/prepsite.f

@@ -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

source_code/version2.2_windows/simul.f

@@ -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

source_code/version2.2_windows/soil_tem.f

@@ -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

source_code/version2.2_windows/sr_forska.f

-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-!       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                                                              

source_code/version2.2_windows/stand_bal.f

@@ -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

source_code/version2.2_windows/win_prepini.f

-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

source_code/version2.2_windows/wpm_wood_proc.f

@@ -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