Line data    Source code

       1             : !> \file mo_integrate.f90
       2             : !> \brief \copybrief mo_integrate
       3             : !> \details \copydetails mo_integrate
       4             : 
       5             : !> \brief Provides integration routines
       6             : !> \details This module provides routine for numerical integration such a Newton-Cotes formulas, etc.
       7             : !> \authors Matthias Cuntz
       8             : !> \date Mar 2013
       9             : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
      10             : !! FORCES is released under the LGPLv3+ license \license_note
      11             : MODULE mo_integrate
      12             : 
      13             :   USE mo_kind, ONLY: i4, sp, dp
      14             : 
      15             :   IMPLICIT NONE
      16             : 
      17             :   PUBLIC :: int_regular ! Integrate regularily spaced data
      18             : 
      19             :   ! ------------------------------------------------------------------
      20             :   !>        \brief Integrate regularily spaced data.
      21             :   !>        \details Integrates regularily spaced data with a 5-point Newton-Cotes formula:
      22             :   !!        \f[ \int y dx = \frac{2}{45} \sum_{i=5,n-4,4} 7 y_{i-4} + 32 y_{i-3} + 12 y_{i-2} + 32 y_{i-1} + 7 y_{i} \f]
      23             :   !!
      24             :   !!        dx=1 if not given.
      25             :   !!
      26             :   !!        \b Example
      27             :   !!        \code{.f90}
      28             :   !!        vec = (/ 1., 2, 3., 4., 5., 6., 7., 8., 9. /)
      29             :   !!        m   = int_regular(vec)
      30             :   !!        -> see also example in test directory
      31             :   !!        \endcode
      32             :   !!
      33             :   !>        \param[in] "real(sp/dp) :: dat(:)"        \f$ y_i \f$ 1D-array with y-values.
      34             :   !>        \param[in] "real(sp/dp) :: dx"             x-spacing (default=1.)
      35             :   !>        \return     real(sp/dp) :: integral — \f$ \int y \f$ integral of y values
      36             : 
      37             :   !>        \author Matthias Cuntz
      38             :   !>        \date Mar 2013
      39             :   INTERFACE int_regular
      40             :      MODULE PROCEDURE int_regular_sp, int_regular_dp
      41             :   END INTERFACE int_regular
      42             : 
      43             :   ! ------------------------------------------------------------------
      44             : 
      45             :   PRIVATE
      46             : 
      47             :   ! ------------------------------------------------------------------
      48             : 
      49             : CONTAINS
      50             : 
      51             :   ! ------------------------------------------------------------------
      52             : 
      53             :   ! integrates tabulated function with equal distant points usign 5-point Newton-Cotes formula
      54          87 :   FUNCTION int_regular_dp(dat, dx)
      55             : 
      56             :     IMPLICIT NONE
      57             : 
      58             :     REAL(dp), DIMENSION(:),           INTENT(IN)  :: dat
      59             :     REAL(dp),               OPTIONAL, INTENT(IN)  :: dx
      60             :     REAL(dp)                                      :: int_regular_dp
      61             : 
      62             :     INTEGER(i4) :: n, n0
      63          87 :     REAL(dp)    :: ddx
      64             : 
      65          87 :     if (size(dat,1) < 5) stop 'Error int_regular_dp: size(dat) < 5'
      66             : 
      67          87 :     if (present(dx)) then
      68          86 :        ddx = dx*2.0_dp/45.0_dp
      69             :     else
      70             :        ddx = 2.0_dp/45.0_dp
      71             :     endif
      72             : 
      73          87 :     n0 = 5
      74          87 :     n  = size(dat,1)
      75             : 
      76          87 :     if (ddx .gt. 0.0_dp) then
      77             :        int_regular_dp = sum( &
      78           0 :             (7.0_dp*(dat(n0-4:n-4:4) + dat(n0:n:4)) + &
      79           0 :             32.0_dp*(dat(n0-3:n-3:4) + dat(n0-1:n-1:4)) + &
      80       33393 :             12.0_dp*dat(n0-2:n-2:4)) )
      81             :        ! to avoid underflow issues
      82          82 :        if ( ddx .lt. 1.0_dp ) then
      83          82 :           if ( int_regular_dp .gt. tiny(1.0_dp)/ddx ) then
      84          82 :              int_regular_dp = ddx * int_regular_dp
      85             :           else
      86             :              int_regular_dp = tiny(1.0_dp)
      87             :           end if
      88             :        else
      89           0 :           int_regular_dp = ddx * int_regular_dp
      90             :        end if
      91             :     else
      92             :        int_regular_dp = 0.0_dp
      93             :     end if
      94             : 
      95          87 :   END FUNCTION int_regular_dp
      96             : 
      97           1 :   FUNCTION int_regular_sp(dat, dx)
      98             : 
      99             :     IMPLICIT NONE
     100             : 
     101             :     REAL(sp), DIMENSION(:),           INTENT(IN)  :: dat
     102             :     REAL(sp),               OPTIONAL, INTENT(IN)  :: dx
     103             :     REAL(sp)                                      :: int_regular_sp
     104             : 
     105             :     INTEGER(i4) :: n, n0
     106           1 :     REAL(sp)    :: ddx
     107             : 
     108           1 :     if (size(dat,1) < 5) stop 'Error int_regular_sp: size(dat) < 5'
     109             : 
     110           1 :     if (present(dx)) then
     111           0 :        ddx = dx*2.0_sp/45.0_sp
     112             :     else
     113             :        ddx = 2.0_sp/45.0_sp
     114             :     endif
     115             : 
     116           1 :     n0 = 5
     117           1 :     n  = size(dat,1)
     118           1 :     if (ddx .gt. 0.0_sp) then
     119             :        int_regular_sp = sum( &
     120           0 :             (7.0_sp*(dat(n0-4:n-4:4) + dat(n0:n:4)) + &
     121           0 :             32.0_sp*(dat(n0-3:n-3:4) + dat(n0-1:n-1:4)) + &
     122           3 :             12.0_sp*dat(n0-2:n-2:4)) )
     123             :        ! to avoid underflow issues
     124           1 :        if ( ddx .lt. 1.0_sp ) then
     125           1 :           if ( int_regular_sp .gt. tiny(1.0_sp)/ddx ) then
     126           1 :              int_regular_sp = ddx * int_regular_sp
     127             :           else
     128             :              int_regular_sp = tiny(1.0_sp)
     129             :           end if
     130             :        else
     131           0 :           int_regular_sp = ddx * int_regular_sp
     132             :        end if
     133             :     else
     134             :        int_regular_sp = 0.0_sp
     135             :     end if
     136             : 
     137          87 :   END FUNCTION int_regular_sp
     138             : 
     139             : END MODULE mo_integrate