!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!       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 für 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                                                              
