7 !! Computation of the positive degree days (PDD) with statistical temperature
8 !! fluctuations; based on semi-analytical solution by Calov and Greve (2005).
9 !! Note that this subroutine uses years as time unit
10 !! (as opposed to the seconds which are usually used by SICOPOLIS).
14 !! Copyright 2009-2013 Reinhard Calov, Ralf Greve
18 !! This file is part of SICOPOLIS.
20 !! SICOPOLIS is free software: you can redistribute it and/or modify
21 !! it under the terms of the GNU General Public License as published by
22 !! the Free Software Foundation, either version 3 of the License, or
23 !! (at your option) any later version.
25 !! SICOPOLIS is distributed in the hope that it will be useful,
26 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
27 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28 !! GNU General Public License for more details.
30 !! You should have received a copy of the GNU General Public License
31 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
33 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
36 !> Computation of the positive degree days (PDD) with statistical temperature
37 !! fluctuations; based on semi-analytical solution by Calov and Greve (2005).
38 !! Note that this subroutine uses years as time unit
39 !! (as opposed to the seconds which are usually used by SICOPOLIS).
40 !<------------------------------------------------------------------------------
41 subroutine pdd(temp_mm, s_stat, ET)
47 real(dp),
dimension(12),
intent(in) :: temp_mm
48 real(dp),
intent(in) :: s_stat
50 real(dp),
intent(out) :: et
53 real(dp) :: inv_sqrt2pi, inv_s_stat, inv_sqrt2
57 real(dp),
parameter :: pi=3.141592653589793_dp
58 real(dp),
parameter :: time_year = 1.0_dp, &
59 time_year_inv = 1.0_dp/time_year, &
60 d_time = 1.0_dp/12.0_dp
62 inv_sqrt2pi = 1.0_dp/sqrt(2.0_dp*pi)
63 inv_s_stat = 1.0_dp/s_stat
64 inv_sqrt2 = 1.0_dp/sqrt(2.0_dp)
70 + ( s_stat*inv_sqrt2pi*exp(-0.5_dp*(temp_mm(n)*inv_s_stat)**2) &
71 + 0.5_dp*temp_mm(n)*
erfcc(-temp_mm(n)*inv_s_stat*inv_sqrt2) ) &
75 et = pdd_sum*time_year_inv