!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*   SR photoper                                                 *!
!*                                                               *!
!*   contains follow global units:                               *!
!* photoper         function for calculation of photoperiod      *!
!* daylength        calculation of day length                    *!
!* avg_sun_incl     Calculates average sun declination for       *!
!                   the season at the given latitude in degrees  *!
!* fixclimscen      subroutine for calculation of delta T and P  *!  
!* glob_rad         Estimation of global radiation from sunshine *!
!* frost_index_total subroutine for calculation of frost index   *!
!*                                                               *!
!*                  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                *!
!*                                                               *!
!*****************************************************************!

REAL FUNCTION PHOTOPER(d,xlatitude)
!     by Thomas Kartschall 8.7.92
!
!     PhotoPeriod -Potential daily Sun Shine Period [h]
!     d           -Ordinal Number of Julian Date [Real!4]
!     latitude    -Latitude by Radiant [Real!4]
!     Northern L>0; Southern L<0
!
!     Polarkreis bei je 66.55° bzw 66°33'36'' N/SL
!
USE data_par
USE data_simul

real d, xlatitude, del, ws, ws2
!
! Equator from 0,2° respectively 0°12'
!
IF (abs(xlatitude).lt.0.0024) then
      photoper=12.0
      return
ENDIF
!
!pole surrounding ab 89,8° bzw 89°48' 
!
IF (xlatitude.ge. 1.567305668)xlatitude= 1.567305668
IF (xlatitude.le.-1.567305668) xlatitude=-1.567305668

      g=2*pi*(d-1.0)/365.25
      del=0.006918-0.399912*cos(g)
      del=del+0.070257*sin(g)-0.006758*cos(g+g)
      del=del+0.000907*sin(g+g)-0.002697*cos(g+g+g)
      del=del+0.00148*sin(g+g+g)
      ws=sin(xlatitude)*sin(del)
      ws2=cos(xlatitude)*cos(del)
!
!polar night duration per day no longer than 24h
!
IF (ws/ws2.ge.1.0) ws=ws2
IF (ws/ws2.le.-1.0) ws=-ws2

      ws=acos(ws/ws2)
      ws=12.*(1.-ws/pi)
!day length is dopple the time between HighNoon and SunRise
      PHOTOPER=2.*(ws)
RETURN
END FUNCTION photoper

!*******************************************************************

FUNCTION DAYLENGTH(doy,hlat)

    USE data_par

    IMPLICIT NONE

    REAL :: hlat
    REAL :: daylength
    INTEGER :: doy

    REAL :: decl,arg

    decl = -23.45*(PI/180)*COS(2.*PI/365.*(doy+10))

    ! latitude is converted to rad
    arg = -TAN(hlat*PI/180.)*TAN(decl);

    IF( arg < -1. ) THEN
        daylength = 24.
    ELSE IF ( arg > 1. ) THEN
        daylength = 0.
    ELSE
        daylength = (24./PI)*ACOS(arg)
    ENDIF

END FUNCTION DAYLENGTH

!*******************************************************************

FUNCTION AVG_SUN_INCL(hlat)

!Calculates average sun declination for the season
! at the given latitude in degrees

  use data_par
  implicit none

  REAL :: avg_sun_incl
  REAL :: hlat, h1, h2, h3
  REAL :: decl, sumbeta, dl, sumdl
  INTEGER :: i, j
  REAL, EXTERNAL :: daylength

  sumdl = 0
  sumbeta = 0

  h1 = sin(PI*hlat/180)
  h2 = cos(PI*hlat/180)

  DO i=120,280,+1
      decl = -23.45*(PI/180)*COS(2.*PI/365.*(i+10))
      dl   = DAYLENGTH(i,hlat)
      ! sun declination at noon
	  h3 = h1*sin(decl)+h2*cos(decl)
	  if(h3.gt.1.) h3 = 1
      avg_sun_incl = 180/PI*asin(h3);

      sumbeta = sumbeta + avg_sun_incl*dl;
      sumdl = sumdl + dl;
  END DO

  avg_sun_incl = sumbeta/sumdl

END FUNCTION AVG_SUN_INCL

!*******************************************************************

SUBROUTINE fixclimscen
! fixclimscen calculates deltaT and deltaPrec for climate change scenarios with
! fixed offsets in temperature and precipitation 
USE data_simul
IMPLICIT NONE
INTEGER   ::  dimTsteps, dimPsteps

! calculations
dimTsteps = 1 + n_T_downsteps + n_T_upsteps
dimPsteps = 1 + n_P_downsteps + n_P_upsteps

deltaT = ((ip-1)/dimPsteps-n_T_downsteps)*step_sum_T
deltaPrec = 1.+((ip-1)-((ip-1)/dimPsteps)*dimPsteps-n_P_downsteps)*step_fac_P

CALL out_scen

END SUBROUTINE fixclimscen

!****************************************************************************

SUBROUTINE glob_rad(sd, iday, xlat, rad)

! Estimation of global radiation from sunshine duration
! (calculation after Angstrom)
implicit none

! input:
integer :: iday        ! actual day
real    :: sd          ! sunshine duration (h)
real    :: xlat        ! latitude

! output:
real    :: rad         ! global radiation  (J/cm2)

! internal variables
real  :: rad_ex  , &   ! extraterrestrical radiation (J/cm2)
         dayl    , &   ! daylength
		 dec     , &   ! declination of sun angle
		 sinld, cosld, tanld, dsinb, dsinbe,  &
         sc, radi, seas
real  :: pi = 3.141592654
real  :: solc = 1367.  ! solar constant (J/(m2*s) 
                       ! after P. Hupfer: "Klimasystem der Erde", 1991                            

!   change of units from degree to radians
pi   = 3.141592654
radi = pi/180. 
 
!   term of seasonality (10 days in front of calendar)
seas = (iday+10.)/365. 
 
!   declination of sun angle
!   (Spitters et al. 1986, equations transformed for use or radians)
dec = -asin(sin(23.45*radi)*cos(2.*pi*seas)) 
 
!   some intermediate values
sinld = sin(xlat*radi)*sin(dec) 
cosld = cos(xlat*radi)*cos(dec) 
tanld = amax1(-1., amin1(1., sinld/cosld)) 
 
!   daylength 
dayl = 12.*(1.+2.*asin(tanld)/pi) 
 
!   integral of sun elevation
dsinb = 3600.*(dayl*sinld+24.*cosld*sqrt(1.-tanld*tanld)/pi) 

!   corrected integral of sun elevation
dsinbe = 3600.*(dayl*(sinld+0.4*(sinld*sinld+cosld*cosld*0.5))  &
         +12.*cosld*(2.+3.*0.4*sinld)*sqrt(1.-tanld*tanld)/pi) 
 
!   intensity of radiation outside the atmosphere
      sc    = solc/(1.-0.016729*cos((360./365.)*(iday-4.)*radi))**2.
      rad_ex = sc*(1.+0.033*cos(2.*pi*iday/365.))*dsinbe
      
!   unit conversion in MJ/m2: rad_ex = rad_ex/1000000. 
!   unit conversion in J/cm2
rad_ex = rad_ex * 0.0001
      
if (sd.ge.0.) then
    rad = (0.231+0.539*sd/dayl)*rad_ex
else   
    write (*, '(A, I3, A)') ' RAD is out of range at day ', iday , &
    ' ,    RAD    will be = 1000 J/cm2!'
end if
           
END	SUBROUTINE glob_rad

!****************************************************************************

subroutine  frost_index_total

use data_frost
use data_simul
use data_stand

implicit none
integer   :: zaehl=0
integer   :: i
integer   :: zaehl1 =0
integer   :: t,m,j
real      :: mean_dnlf
real      :: mean_tminmay
integer   :: mean_date_lf
integer   :: mean_date_lftot
real      :: mean_dnlf_sp
real      :: mean_tminmay_sp
integer   :: mean_date_lf_sp
real      :: mean_anzdlf
real      :: mean_sumtlf
integer   :: ind1, ind2, ind3, ind4, ind5
integer   :: ind1_sp

zaehl=0
mean_tminmay = 0.
mean_date_lf = 0
mean_date_lftot = 0
mean_dnlf = 0
mean_dnlf_sp = 0
mean_anzdlf = 0
mean_sumtlf = 0

do i =1,year
 if(tminmay_ann(i).ne.0) then
   zaehl = zaehl +1
   mean_tminmay= mean_tminmay+tminmay_ann(i)
 end if
end do  

if(zaehl.ne.0) then
   mean_tminmay = mean_tminmay/zaehl
else
  mean_tminmay = 0.
end if

do i=1,year
   mean_anzdlf = mean_anzdlf + anzdlf(i)
   mean_sumtlf = mean_sumtlf + sumtlf(i)
end do

mean_anzdlf = mean_anzdlf/year
mean_sumtlf = mean_sumtlf/year

zaehl=0
do i =1,year
 if(date_lftot(i).ne.0) then
   zaehl = zaehl +1
   mean_date_lftot = mean_date_lftot + date_lftot(i)
 end if
end do  

if(zaehl.ne.0) then
   mean_date_lftot = mean_date_lftot/zaehl
else
  mean_date_lftot = 0.
end if
mean_dnlf = 0.
zaehl=0
do i =1,year
   if(dnlf(i).ne.0) then
     mean_dnlf = mean_dnlf + dnlf(i)
     zaehl = zaehl +1
   end if
end do
if(zaehl.ne.0) then
   mean_dnlf = mean_dnlf/zaehl
else
   mean_dnlf = 0
endif

zaehl=0
do i =1,year
   if(date_lf(i).ne.0) then
       mean_date_lf = mean_date_lf + date_lf(i)
       zaehl = zaehl +1
   end if
enddo

if(zaehl.ne.0) then
  mean_date_lf = mean_date_lf/zaehl
else
 mean_date_lf = 0
end if

zaehl1=0
do i =1,year
 if(dnlf_sp(i).ne.0) then
   zaehl1 = zaehl1 +1
   mean_dnlf_sp = mean_dnlf_sp + dnlf_sp(i)
 end if
enddo

if (zaehl1.ne.0) then
   mean_dnlf_sp = mean_dnlf_sp/zaehl1
 else
   mean_dnlf_sp = 0
endif

if (mean_dnlf.le.2.5 .and. mean_tminmay.ge. -1.5 .and.tminmay.ge.-5.0 .and. mean_date_lf.lt.130 .and. dlfabs .lt. 156) lfind=1
if (mean_dnlf.ge.2.6 .and. mean_dnlf .le.3.5 .and. mean_tminmay.ge. -2.0 .and. mean_tminmay.lt.-1.5 .and. tminmay .ge.-6. .and. mean_date_lf .lt.135 .and. dlfabs .lt.161) lfind=2
if (mean_dnlf.gt.3.5 .and. mean_dnlf .le.4.5 .and. mean_tminmay.ge. -2.5 .and. mean_tminmay.lt.-2.0 .and. tminmay .ge.-6. .and. mean_date_lf .ge.135 .and. mean_date_lf .le. 140 .and. dlfabs .ge.162 .and. dlfabs.le.166) lfind=3
if (mean_dnlf.gt.4.5 .and. mean_dnlf .le.5.0 .and. mean_tminmay.ge. -3.0 .and. mean_tminmay.lt.-2.5 .and. tminmay .ge.-7. .and. mean_date_lf .ge.141 .and. mean_date_lf .le. 145 .and. dlfabs .ge.167 .and. dlfabs.le.171) lfind=4
if (mean_dnlf.gt.5.10 .and. mean_dnlf .le.5.5 .and. mean_tminmay.ge. -3.5 .and. mean_tminmay.lt.-3.0 .and. tminmay .ge.-8. .and. mean_date_lf .ge.141 .and. mean_date_lf .le. 145 .and. dlfabs .ge.172 .and. dlfabs.le.176) lfind=5
if (mean_dnlf.gt.5.5  .and. mean_tminmay.lt.-3.5 .and. tminmay .le.-8. .and. mean_date_lf .gt.145 .and. dlfabs .gt.176) lfind=6

! index of number of late frost days since beginning of vegetation period
if (mean_dnlf.le.2.5) then
   ind1 = 1
else if(mean_dnlf.le.3.5) then
   ind1 = 2
else if (mean_dnlf.le.4.5) then
   ind1 = 3
else if (mean_dnlf.le.5.0) then
   ind1 = 4
else if (mean_dnlf.le.5.5) then
   ind1 = 5
else
   ind1 = 6
endif
! index of number of late frost days since beginning of bud burst
if (mean_dnlf_sp .le. 2.5) then   
   ind1_sp= 1
else if(mean_dnlf_sp.le.3.5) then
   ind1_sp = 2
else if (mean_dnlf_sp.le.4.5) then
   ind1_sp = 3
else if (mean_dnlf.le.5.0) then
   ind1_sp = 4
else if (mean_dnlf_sp.le.5.5) then
   ind1_sp = 5
else
   ind1_sp = 6
endif
! index of mean minimum may temperature
if(mean_tminmay.ge. -1.5) then
  ind2 = 1
else if (mean_tminmay.ge. -2.0) then
  ind2 = 2
else if (mean_tminmay.ge. -2.5) then
 ind2 = 3
else if (mean_tminmay.ge. -3.0) then
 ind2 = 4
else if (mean_tminmay.ge. -3.5) then
 ind2 = 5
else
 ind2 =6
endif
! index of absolute minimum may temperature
if(tminmay.ge.-5.0) then
  ind3 = 1
else if(tminmay.ge.-6.0 .and. ind2 .le.2) then
  ind3 = 2
else if (tminmay.ge.-6.0 .and. ind2 .le.3) then
  ind3 =3
else if (tminmay.ge.-7.0) then
  ind3 = 4
else if (tminmay.ge.-8.0) then
  ind3 = 5
else
  ind3 = 6
end if
! index of mean date(number of the year) of late frost
if (mean_date_lf.lt.130) then
  ind4 = 1
else if (mean_date_lf.lt.135) then
  ind4 = 2
else if (mean_date_lf.le.140 ) then
  ind4 = 3
else if (mean_date_lf.le.145 .and. ind2.le.4) then
  ind4 = 4
else if(mean_date_lf.le.145 .and. ind2.le.5) then
  ind4 = 5
  else
   ind4 = 6
endif
! absolute last late frost (numbedr of the year)
if (dlfabs .lt. 156) then
  ind5 = 1
else if (dlfabs .lt. 161) then
  ind5 = 2
else if (dlfabs .le. 162) then
  ind5 =3
else if (dlfabs .le. 171) then
  ind5 = 4
else if (dlfabs .le. 176) then
  ind5 = 5
else
  ind5 =6
endif  

mlfind = real((ind1 + ind2 + ind3 + ind4 + ind5)/5)
mlfind_sp = (ind1_sp + ind2 + ind3 + ind4 + ind5)/5
  
  if(waldtyp.eq. 10 .or. waldtyp .eq. 40 .or. waldtyp .eq.90) mlfind_sp = 0
      
end subroutine frost_index_total






