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