README.md

@@ -18,8 +18,10 @@ collaborative work, see the AUTHOR file for details.
     * manuals: about "how to start" and " generate exe"  
     * descriptions/additional: background information documents 
     * descriptions/Fact_sheets: short summary of the model in English and German
+    * source_code/additional: overview on variables with their names in the 4C description and the names in the source code
+                              flow chart of the model 4C 
     * test_data_sets: test the correct function mode of the model with complete sets of input data and control files 
-    * publications: all relevant publications regarding the model 4C
+    * publications: relevant publications regarding the model 4C
 
 The source code is distributed via a git repository at: 
 https://gitlab.pik-potsdam.de/foresee/4C and is ready to use.
@@ -27,7 +29,7 @@ The code language is Fortran90. Since the model development is
 finished, there is no further development and no further support in model 
 download, setup, development or application. No merge or other development 
 requests will be handeled. Therefore there are no code style recommendations or
-instructions on how to create source code ammendments.
+instructions on how to create source code amendments.
 
 
 

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