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