pdd.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  p d d
!
!> @file
!!
!! Computation of the positive degree days (PDD) with statistical temperature
!! fluctuations; based on semi-analytical solution by Calov and Greve (2005).
!! Note that this subroutine uses years as time unit
!! (as opposed to the seconds which are usually used by SICOPOLIS).
!!
!! @section Copyright
!!
!! Copyright 2009-2013 Reinhard Calov, Ralf Greve
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Computation of the positive degree days (PDD) with statistical temperature
!! fluctuations; based on semi-analytical solution by Calov and Greve (2005).
!! Note that this subroutine uses years as time unit
!! (as opposed to the seconds which are usually used by SICOPOLIS).
!<------------------------------------------------------------------------------
subroutine pdd(temp_mm, s_stat, ET)

use sico_types

implicit none

real(dp), dimension(12), intent(in) :: temp_mm
real(dp), intent(in) :: s_stat

real(dp), intent(out) :: et

integer(i4b) :: n
real(dp) :: inv_sqrt2pi, inv_s_stat, inv_sqrt2
real(dp) :: pdd_sum
real(dp) :: erfcc

real(dp), parameter :: pi=3.141592653589793_dp
real(dp), parameter :: time_year     = 1.0_dp, &        ! period 1 year
                       time_year_inv = 1.0_dp/time_year, &
                       d_time        = 1.0_dp/12.0_dp   ! time-step 1 month

inv_sqrt2pi = 1.0_dp/sqrt(2.0_dp*pi)
inv_s_stat  = 1.0_dp/s_stat
inv_sqrt2   = 1.0_dp/sqrt(2.0_dp)

pdd_sum = 0.0_dp

do n=1, 12   ! month counter
   pdd_sum = pdd_sum &
             + ( s_stat*inv_sqrt2pi*exp(-0.5_dp*(temp_mm(n)*inv_s_stat)**2) &
                 + 0.5_dp*temp_mm(n)*erfcc(-temp_mm(n)*inv_s_stat*inv_sqrt2) ) &
               *d_time   ! positive degree days (in a * deg C)
end do

et   = pdd_sum*time_year_inv   ! temperature excess   (in deg C)

end subroutine pdd
!